!>\file funcphys.f90 
!! This file includes API for basic thermodynamic physics.

!>\defgroup func_phys GFS Physics Function Module
!! This module provides API for computing basic thermodynamic physics
!! functions. 

!> This module provides an Application Program Interface (API) for computing
!! basic thermodynamic physics functions, in particular:
!! -# saturation vapor pressure as a function of temperature;
!! -# dewpoint temperature as a function of vapor pressure;
!! -# equivalent potential temperature as a function of temperature and 
!! scaled pressure to the kappa power;
!! -# temperature and specific humidity along a moist adiabat as functions
!! of equivalent potential temperature and scaled pressure to the kappa power;
!! -# scaled pressure to the kappa power as a function of pressure, and
!! -# temperature at the lifting condensation level as a function of temperature
!! and dewpoint depression.
!!
!! The entry points required to set up lookup tables start with a "g".
!! All the other entry points are functions starting with an "f" or 
!! are subroutines starting with an "s". These other functions and subroutines
!! are elemental; that is, they return a scalar if they are passed only scalars,
!! but they return an array if they are passed an array. These other functions
!! and subroutines can in inlined, too.
module funcphys
!$$$  Module Documentation Block
!
! Module:    funcphys        API for basic thermodynamic physics
!   Author: Iredell          Org: W/NX23     Date: 1999-03-01
!
! Abstract: This module provides an Application Program Interface
!   for computing basic thermodynamic physics functions, in particular
!     (1) saturation vapor pressure as a function of temperature,
!     (2) dewpoint temperature as a function of vapor pressure,
!     (3) equivalent potential temperature as a function of temperature
!         and scaled pressure to the kappa power,
!     (4) temperature and specific humidity along a moist adiabat
!         as functions of equivalent potential temperature and
!         scaled pressure to the kappa power,
!     (5) scaled pressure to the kappa power as a function of pressure, and
!     (6) temperature at the lifting condensation level as a function
!         of temperature and dewpoint depression.
!   The entry points required to set up lookup tables start with a "g".
!   All the other entry points are functions starting with an "f" or
!   are subroutines starting with an "s".  These other functions and
!   subroutines are elemental; that is, they return a scalar if they
!   are passed only scalars, but they return an array if they are passed
!   an array.  These other functions and subroutines can be inlined, too.
!
! Program History Log:
!   1999-03-01  Mark Iredell
!   1999-10-15  Mark Iredell  SI unit for pressure (Pascals)
!   2001-02-26  Mark Iredell  Ice phase changes of Hong and Moorthi
!
! Public Variables:
!   krealfp         Integer parameter kind or length of reals (=kind_phys)
!
! Public Subprograms:
!   gpvsl            Compute saturation vapor pressure over liquid table
!
!   fpvsl           Elementally compute saturation vapor pressure over liquid
!     function result Real(krealfp) saturation vapor pressure in Pascals
!     t               Real(krealfp) temperature in Kelvin
!
!   fpvslq          Elementally compute saturation vapor pressure over liquid
!     function result Real(krealfp) saturation vapor pressure in Pascals
!     t               Real(krealfp) temperature in Kelvin
!
!   fpvslx          Elementally compute saturation vapor pressure over liquid
!     function result Real(krealfp) saturation vapor pressure in Pascals
!     t               Real(krealfp) temperature in Kelvin
!
!   gpvsi            Compute saturation vapor pressure over ice table
!
!   fpvsi           Elementally compute saturation vapor pressure over ice
!     function result Real(krealfp) saturation vapor pressure in Pascals
!     t               Real(krealfp) temperature in Kelvin
!
!   fpvsiq          Elementally compute saturation vapor pressure over ice
!     function result Real(krealfp) saturation vapor pressure in Pascals
!     t               Real(krealfp) temperature in Kelvin
!
!   fpvsix          Elementally compute saturation vapor pressure over ice
!     function result Real(krealfp) saturation vapor pressure in Pascals
!     t               Real(krealfp) temperature in Kelvin
!
!   gpvs            Compute saturation vapor pressure table
!
!   fpvs            Elementally compute saturation vapor pressure
!     function result Real(krealfp) saturation vapor pressure in Pascals
!     t               Real(krealfp) temperature in Kelvin
!
!   fpvsq           Elementally compute saturation vapor pressure
!     function result Real(krealfp) saturation vapor pressure in Pascals
!     t               Real(krealfp) temperature in Kelvin
!
!   fpvsx           Elementally compute saturation vapor pressure
!     function result Real(krealfp) saturation vapor pressure in Pascals
!     t               Real(krealfp) temperature in Kelvin
!
!   gtdpl           Compute dewpoint temperature over liquid table
!
!   ftdpl           Elementally compute dewpoint temperature over liquid
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   ftdplq          Elementally compute dewpoint temperature over liquid
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   ftdplx          Elementally compute dewpoint temperature over liquid
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   ftdplxg         Elementally compute dewpoint temperature over liquid
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     t               Real(krealfp) guess dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   gtdpi           Compute dewpoint temperature table over ice
!
!   ftdpi           Elementally compute dewpoint temperature over ice
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   ftdpiq          Elementally compute dewpoint temperature over ice
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   ftdpix          Elementally compute dewpoint temperature over ice
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   ftdpixg         Elementally compute dewpoint temperature over ice
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     t               Real(krealfp) guess dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   gtdp            Compute dewpoint temperature table
!
!   ftdp            Elementally compute dewpoint temperature
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   ftdpq           Elementally compute dewpoint temperature
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   ftdpx           Elementally compute dewpoint temperature
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   ftdpxg          Elementally compute dewpoint temperature
!     function result Real(krealfp) dewpoint temperature in Kelvin
!     t               Real(krealfp) guess dewpoint temperature in Kelvin
!     pv              Real(krealfp) vapor pressure in Pascals
!
!   gthe            Compute equivalent potential temperature table
!
!   fthe            Elementally compute equivalent potential temperature
!     function result Real(krealfp) equivalent potential temperature in Kelvin
!     t               Real(krealfp) LCL temperature in Kelvin
!     pk              Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
!
!   ftheq           Elementally compute equivalent potential temperature
!     function result Real(krealfp) equivalent potential temperature in Kelvin
!     t               Real(krealfp) LCL temperature in Kelvin
!     pk              Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
!
!   fthex           Elementally compute equivalent potential temperature
!     function result Real(krealfp) equivalent potential temperature in Kelvin
!     t               Real(krealfp) LCL temperature in Kelvin
!     pk              Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
!
!   gtma            Compute moist adiabat tables
!
!   stma            Elementally compute moist adiabat temperature and moisture
!     the             Real(krealfp) equivalent potential temperature in Kelvin
!     pk              Real(krealfp) pressure over 1e5 Pa to the kappa power
!     tma             Real(krealfp) parcel temperature in Kelvin
!     qma             Real(krealfp) parcel specific humidity in kg/kg
!
!   stmaq           Elementally compute moist adiabat temperature and moisture
!     the             Real(krealfp) equivalent potential temperature in Kelvin
!     pk              Real(krealfp) pressure over 1e5 Pa to the kappa power
!     tma             Real(krealfp) parcel temperature in Kelvin
!     qma             Real(krealfp) parcel specific humidity in kg/kg
!
!   stmax           Elementally compute moist adiabat temperature and moisture
!     the             Real(krealfp) equivalent potential temperature in Kelvin
!     pk              Real(krealfp) pressure over 1e5 Pa to the kappa power
!     tma             Real(krealfp) parcel temperature in Kelvin
!     qma             Real(krealfp) parcel specific humidity in kg/kg
!
!   stmaxg          Elementally compute moist adiabat temperature and moisture
!     tg              Real(krealfp) guess parcel temperature in Kelvin
!     the             Real(krealfp) equivalent potential temperature in Kelvin
!     pk              Real(krealfp) pressure over 1e5 Pa to the kappa power
!     tma             Real(krealfp) parcel temperature in Kelvin
!     qma             Real(krealfp) parcel specific humidity in kg/kg
!
!   gpkap           Compute pressure to the kappa table
!
!   fpkap           Elementally raise pressure to the kappa power.
!     function result Real(krealfp) p over 1e5 Pa to the kappa power
!     p               Real(krealfp) pressure in Pascals
!
!   fpkapq          Elementally raise pressure to the kappa power.
!     function result Real(krealfp) p over 1e5 Pa to the kappa power
!     p               Real(krealfp) pressure in Pascals
!
!   fpkapo          Elementally raise pressure to the kappa power.
!     function result Real(krealfp) p over 1e5 Pa to the kappa power
!     p               Real(krealfp) surface pressure in Pascals
!
!   fpkapx          Elementally raise pressure to the kappa power.
!     function result Real(krealfp) p over 1e5 Pa to the kappa power
!     p               Real(krealfp) pressure in Pascals
!
!   grkap           Compute pressure to the 1/kappa table
!
!   frkap           Elementally raise pressure to the 1/kappa power.
!     function result Real(krealfp) pressure in Pascals
!     pkap            Real(krealfp) p over 1e5 Pa to the 1/kappa power
!
!   frkapq          Elementally raise pressure to the kappa power.
!     function result Real(krealfp) pressure in Pascals
!     pkap            Real(krealfp) p over 1e5 Pa to the kappa power
!
!   frkapx          Elementally raise pressure to the kappa power.
!     function result Real(krealfp) pressure in Pascals
!     pkap            Real(krealfp) p over 1e5 Pa to the kappa power
!
!   gtlcl           Compute LCL temperature table
!
!   ftlcl           Elementally compute LCL temperature.
!     function result Real(krealfp) temperature at the LCL in Kelvin
!     t               Real(krealfp) temperature in Kelvin
!     tdpd            Real(krealfp) dewpoint depression in Kelvin
!
!   ftlclq          Elementally compute LCL temperature.
!     function result Real(krealfp) temperature at the LCL in Kelvin
!     t               Real(krealfp) temperature in Kelvin
!     tdpd            Real(krealfp) dewpoint depression in Kelvin
!
!   ftlclo          Elementally compute LCL temperature.
!     function result Real(krealfp) temperature at the LCL in Kelvin
!     t               Real(krealfp) temperature in Kelvin
!     tdpd            Real(krealfp) dewpoint depression in Kelvin
!
!   ftlclx          Elementally compute LCL temperature.
!     function result Real(krealfp) temperature at the LCL in Kelvin
!     t               Real(krealfp) temperature in Kelvin
!     tdpd            Real(krealfp) dewpoint depression in Kelvin
!
!   gfuncphys       Compute all physics function tables
!
! Attributes:
!   Language: Fortran 90
!
!$$$
  use machine,only:kind_phys,r8=>kind_dbl_prec,r4=>kind_sngl_prec
  use physcons
  implicit none
  private
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Public Variables
! integer,public,parameter:: krealfp=selected_real_kind(15,45)
  integer,public,parameter:: krealfp=kind_phys          !< Integer parameter kind or length of reals
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Private Variables
  real(krealfp),parameter:: psatb=con_psat*1.e-5
  integer,parameter:: nxpvsl=7501
  real(krealfp) c1xpvsl,c2xpvsl,tbpvsl(nxpvsl)
  integer,parameter:: nxpvsi=7501
  real(krealfp) c1xpvsi,c2xpvsi,tbpvsi(nxpvsi)
  integer,parameter:: nxpvs=7501
  real(krealfp) c1xpvs,c2xpvs,tbpvs(nxpvs)
  integer,parameter:: nxtdpl=5001
  real(krealfp) c1xtdpl,c2xtdpl,tbtdpl(nxtdpl)
  integer,parameter:: nxtdpi=5001
  real(krealfp) c1xtdpi,c2xtdpi,tbtdpi(nxtdpi)
  integer,parameter:: nxtdp=5001
  real(krealfp) c1xtdp,c2xtdp,tbtdp(nxtdp)
  integer,parameter:: nxthe=241,nythe=151
  real(krealfp) c1xthe,c2xthe,c1ythe,c2ythe,tbthe(nxthe,nythe)
  integer,parameter:: nxma=151,nyma=121
  real(krealfp) c1xma,c2xma,c1yma,c2yma,tbtma(nxma,nyma),tbqma(nxma,nyma)
! integer,parameter:: nxpkap=5501
  integer,parameter:: nxpkap=11001
  real(krealfp) c1xpkap,c2xpkap,tbpkap(nxpkap)
  integer,parameter:: nxrkap=5501
  real(krealfp) c1xrkap,c2xrkap,tbrkap(nxrkap)
  integer,parameter:: nxtlcl=151,nytlcl=61
  real(krealfp) c1xtlcl,c2xtlcl,c1ytlcl,c2ytlcl,tbtlcl(nxtlcl,nytlcl)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Public Subprograms
  public gpvsl,fpvsl,fpvslq,fpvslx
  public gpvsi,fpvsi,fpvsiq,fpvsix
  public gpvs,fpvs,fpvsq,fpvsx
  public gtdpl,ftdpl,ftdplq,ftdplx,ftdplxg
  public gtdpi,ftdpi,ftdpiq,ftdpix,ftdpixg
  public gtdp,ftdp,ftdpq,ftdpx,ftdpxg
  public gthe,fthe,ftheq,fthex
  public gtma,stma,stmaq,stmax,stmaxg
  public gpkap,fpkap,fpkapq,fpkapo,fpkapx
  public grkap,frkap,frkapq,frkapx
  public gtlcl,ftlcl,ftlclq,ftlclo,ftlclx
  public gfuncphys

  interface fpvsl
     module procedure fpvsl_r4, fpvsl_r8 
  end interface fpvsl
  interface fpvsi
     module procedure fpvsi_r4, fpvsi_r8 
  end interface fpvsi
contains
!-------------------------------------------------------------------------------
!> This subroutine computes saturation vapor pressure table as a function of
!! temperature for the table lookup function fpval. Exact saturation vapor
!! pressures are calculated in subprogram fpvslx(). The current implementation
!! computes a table with a length of 7501 for temperature ranging from 180. to
!! 330. Kelvin.
  subroutine gpvsl
!$$$     Subprogram Documentation Block
!
! Subprogram: gpvsl        Compute saturation vapor pressure table over liquid
!   Author: N Phillips            W/NMC2X2   Date: 30 dec 82
!
! Abstract: Computes saturation vapor pressure table as a function of
!   temperature for the table lookup function fpvsl.
!   Exact saturation vapor pressures are calculated in subprogram fpvslx.
!   The current implementation computes a table with a length
!   of 7501 for temperatures ranging from 180. to 330. Kelvin.
!
! Program History Log:
!   91-05-07  Iredell
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
!
! Usage:  call gpvsl
!
! Subprograms called:
!   (fpvslx)   inlinable function to compute saturation vapor pressure
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx
    real(krealfp) xmin,xmax,xinc,x,t
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=180.0_krealfp
    xmax=330.0_krealfp
    xinc=(xmax-xmin)/(nxpvsl-1)
!   c1xpvsl=1.-xmin/xinc
    c2xpvsl=1./xinc
    c1xpvsl=1.-xmin*c2xpvsl
    do jx=1,nxpvsl
      x=xmin+(jx-1)*xinc
      t=x
      tbpvsl(jx)=fpvslx(t)
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This funtion computes saturation vapor pressure from the temperature.
!! A linear interpolation is done between values in a lookup table computed
!! in gpvsl(). See documentation for fpvslx() for details. Input values
!! outside table range are reset to table extrema. 
!>\author N phillips

  elemental function fpvsl_r4(t)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvsl        Compute saturation vapor pressure over liquid
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute saturation vapor pressure from the temperature.
!   A linear interpolation is done between values in a lookup table
!   computed in gpvsl. See documentation for fpvslx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is almost 6 decimal places.
!   On the Cray, fpvsl is about 4 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
!
! Usage:   pvsl=fpvsl(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvsl      Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(r4) fpvsl_r4
    real(r4),intent(in):: t
    integer jx
    real(r4) xj
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xpvsl+c2xpvsl*t,1._r4),real(nxpvsl,r4))
    jx=min(xj,nxpvsl-1._r4)
    fpvsl_r4=tbpvsl(jx)+(xj-jx)*(tbpvsl(jx+1)-tbpvsl(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function fpvsl_r4

  elemental function fpvsl_r8(t)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvsl        Compute saturation vapor pressure over liquid
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute saturation vapor pressure from the temperature.
!   A linear interpolation is done between values in a lookup table
!   computed in gpvsl. See documentation for fpvslx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is almost 6 decimal places.
!   On the Cray, fpvsl is about 4 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
!
! Usage:   pvsl=fpvsl(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvsl      Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(r8) fpvsl_r8
    real(r8),intent(in):: t
    integer jx
    real(r8) xj
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xpvsl+c2xpvsl*t,1._r8),real(nxpvsl,r8))
    jx=min(xj,nxpvsl-1._r8)
    fpvsl_r8=tbpvsl(jx)+(xj-jx)*(tbpvsl(jx+1)-tbpvsl(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function fpvsl_r8



!-------------------------------------------------------------------------------
!> This function computes saturation vapor pressure from the temperature. 
!! A quadratic interpolation is done between values in a lookup table 
!! computed in gpvsl(). See documentaion for fpvslx() for details.
!! Input values outside table range are reset to table extrema.
  elemental function fpvslq(t)
!>\author N Phillips
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvslq       Compute saturation vapor pressure over liquid
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute saturation vapor pressure from the temperature.
!   A quadratic interpolation is done between values in a lookup table
!   computed in gpvsl. See documentation for fpvslx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is almost 9 decimal places.
!   On the Cray, fpvslq is about 3 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             quadratic interpolation
! 1999-03-01  Iredell             f90 module
!
! Usage:   pvsl=fpvslq(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvslq     Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpvslq
    real(krealfp),intent(in):: t
    integer jx
    real(krealfp) xj,dxj,fj1,fj2,fj3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xpvsl+c2xpvsl*t,1._krealfp),real(nxpvsl,krealfp))
    jx=min(max(nint(xj),2),nxpvsl-1)
    dxj=xj-jx
    fj1=tbpvsl(jx-1)
    fj2=tbpvsl(jx)
    fj3=tbpvsl(jx+1)
    fpvslq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function exactly computes saturation vapor pressure from temperature.
!! The water model assumes a perfect gas, constant specific heats
!! for gas and liquid, and neglects the volume of the liquid.
!! The model does account for the variation of the latent heat 
!! of condensation with temperature. The ice option is not included.
!! The Clausius-Clapeyron equation is integrated from the triple point
!! to get the formula:
!!\n pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
!!\n where tr is ttp/t and other values are physical constants.
!! This function should be expanded inline in the calling routine.
!>\author N Phillips
!>\param[in] t       real, temperature in Kelvin
!\param[out] fpvslx real, saturation vapor pressure in Pascals
  elemental function fpvslx(t)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvslx       Compute saturation vapor pressure over liquid
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Exactly compute saturation vapor pressure from temperature.
!   The water model assumes a perfect gas, constant specific heats
!   for gas and liquid, and neglects the volume of the liquid.
!   The model does account for the variation of the latent heat
!   of condensation with temperature.  The ice option is not included.
!   The Clausius-Clapeyron equation is integrated from the triple point
!   to get the formula
!       pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
!   where tr is ttp/t and other values are physical constants.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
!
! Usage:   pvsl=fpvslx(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvslx     Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpvslx
    real(krealfp),intent(in):: t
    real(krealfp),parameter:: dldt=con_cvap-con_cliq
    real(krealfp),parameter:: heat=con_hvap
    real(krealfp),parameter:: xpona=-dldt/con_rv
    real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
    real(krealfp) tr
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    tr=con_ttp/t
    fpvslx=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This subroutine computes saturation vapor pressure table as a function of 
!! temperature for the table lookup function fpvsi(). Exact saturation vapor
!! pressures are calculated in subprogram fpvsix(). The current implementation
!! computes a table with a length of 7501 for temperatures ranging from 180. 
!! to 330. Kelvin.
!>\author N Phillips
  subroutine gpvsi
!$$$     Subprogram Documentation Block
!
! Subprogram: gpvsi        Compute saturation vapor pressure table over ice
!   Author: N Phillips            W/NMC2X2   Date: 30 dec 82
!
! Abstract: Computes saturation vapor pressure table as a function of
!   temperature for the table lookup function fpvsi.
!   Exact saturation vapor pressures are calculated in subprogram fpvsix.
!   The current implementation computes a table with a length
!   of 7501 for temperatures ranging from 180. to 330. Kelvin.
!
! Program History Log:
!   91-05-07  Iredell
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:  call gpvsi
!
! Subprograms called:
!   (fpvsix)   inlinable function to compute saturation vapor pressure
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx
    real(krealfp) xmin,xmax,xinc,x,t
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=180.0_krealfp
    xmax=330.0_krealfp
    xinc=(xmax-xmin)/(nxpvsi-1)
!   c1xpvsi=1.-xmin/xinc
    c2xpvsi=1./xinc
    c1xpvsi=1.-xmin*c2xpvsi
    do jx=1,nxpvsi
      x=xmin+(jx-1)*xinc
      t=x
      tbpvsi(jx)=fpvsix(t)
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This function computes saturation vapor pressure from the temperature.
!! A linear interpolation is done between values in a lookup table 
!! computed in gpvsi(). See documentation for fpvsix() for details.
!! Input values outside table range are reset to table extrema.
!>\author N Phillips

  elemental function fpvsi_r4(t)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvsi        Compute saturation vapor pressure over ice
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute saturation vapor pressure from the temperature.
!   A linear interpolation is done between values in a lookup table
!   computed in gpvsi. See documentation for fpvsix for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is almost 6 decimal places.
!   On the Cray, fpvsi is about 4 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   pvsi=fpvsi(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvsi      Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(r4) fpvsi_r4
    real(r4),intent(in):: t
    integer jx
    real(r4) xj
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xpvsi+c2xpvsi*t,1._r4),real(nxpvsi,r4))
    jx=min(xj,nxpvsi-1._r4)
    fpvsi_r4=tbpvsi(jx)+(xj-jx)*(tbpvsi(jx+1)-tbpvsi(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function fpvsi_r4

  elemental function fpvsi_r8(t)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvsi        Compute saturation vapor pressure over ice
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute saturation vapor pressure from the temperature.
!   A linear interpolation is done between values in a lookup table
!   computed in gpvsi. See documentation for fpvsix for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is almost 6 decimal places.
!   On the Cray, fpvsi is about 4 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   pvsi=fpvsi(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvsi      Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(r8) fpvsi_r8
    real(r8),intent(in):: t
    integer jx
    real(r8) xj
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xpvsi+c2xpvsi*t,1._r8),real(nxpvsi,r8))
    jx=min(xj,nxpvsi-1._r8)
    fpvsi_r8=tbpvsi(jx)+(xj-jx)*(tbpvsi(jx+1)-tbpvsi(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function fpvsi_r8

!-------------------------------------------------------------------------------
!> This function computes saturation vapor pressure from the temperature.
!! A quadratic interpolation is done between values in a lookup table
!! computed in gpvsi(). See documentation for fpvsix() for details.
!! Input values outside table range are reset to table extrema.
!>\author N Phillips
  elemental function fpvsiq(t)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvsiq       Compute saturation vapor pressure over ice
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute saturation vapor pressure from the temperature.
!   A quadratic interpolation is done between values in a lookup table
!   computed in gpvsi. See documentation for fpvsix for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is almost 9 decimal places.
!   On the Cray, fpvsiq is about 3 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             quadratic interpolation
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   pvsi=fpvsiq(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvsiq     Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpvsiq
    real(krealfp),intent(in):: t
    integer jx
    real(krealfp) xj,dxj,fj1,fj2,fj3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xpvsi+c2xpvsi*t,1._krealfp),real(nxpvsi,krealfp))
    jx=min(max(nint(xj),2),nxpvsi-1)
    dxj=xj-jx
    fj1=tbpvsi(jx-1)
    fj2=tbpvsi(jx)
    fj3=tbpvsi(jx+1)
    fpvsiq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This funtion exactly computes saturation vapor pressure from temperature.
!! The water model assumes a perfect gas, constant specific heats
!! for gas and ice, and neglects the volume of the ice. The model does
!! account for the variation of the latent heat of condensation with temperature.
!! The liquid option is not included. The Clausius- Clapeyron equation is 
!! integrated from the triple point to get the formula:
!!\n pvsi=con_psat*(tr**xa)*exp(xb*(1.-tr))
!!\n where tr is ttp/t and other values are physical constants. 
!! This function should be expanded inline in the calling routine.
!>\param[in]  t        real, temperature in Kelvin
!\param[out] fpvsix   real, saturation vapor pressure in Pascals
  elemental function fpvsix(t)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvsix       Compute saturation vapor pressure over ice
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Exactly compute saturation vapor pressure from temperature.
!   The water model assumes a perfect gas, constant specific heats
!   for gas and ice, and neglects the volume of the ice.
!   The model does account for the variation of the latent heat
!   of condensation with temperature.  The liquid option is not included.
!   The Clausius-Clapeyron equation is integrated from the triple point
!   to get the formula
!       pvsi=con_psat*(tr**xa)*exp(xb*(1.-tr))
!   where tr is ttp/t and other values are physical constants.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   pvsi=fpvsix(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvsix     Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpvsix
    real(krealfp),intent(in):: t
    real(krealfp),parameter:: dldt=con_cvap-con_csol
    real(krealfp),parameter:: heat=con_hvap+con_hfus
    real(krealfp),parameter:: xpona=-dldt/con_rv
    real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
    real(krealfp) tr
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    tr=con_ttp/t
    fpvsix=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This subroutine computes saturation vapor pressure table as a function of
!! temperature for the table lookup function fpvs().
!! Exact saturation vapor pressures are calculated in subprogram fpvsx().
!! The current implementation computes a table with a length
!! of 7501 for temperatures ranging from 180. to 330. Kelvin. 
  subroutine gpvs
!$$$     Subprogram Documentation Block
!
! Subprogram: gpvs         Compute saturation vapor pressure table
!   Author: N Phillips            W/NMC2X2   Date: 30 dec 82
!
! Abstract: Computes saturation vapor pressure table as a function of
!   temperature for the table lookup function fpvs.
!   Exact saturation vapor pressures are calculated in subprogram fpvsx.
!   The current implementation computes a table with a length
!   of 7501 for temperatures ranging from 180. to 330. Kelvin.
!
! Program History Log:
!   91-05-07  Iredell
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:  call gpvs
!
! Subprograms called:
!   (fpvsx)    inlinable function to compute saturation vapor pressure
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx
    real(krealfp) xmin,xmax,xinc,x,t
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=180.0_krealfp
    xmax=330.0_krealfp
    xinc=(xmax-xmin)/(nxpvs-1)
!   c1xpvs=1.-xmin/xinc
    c2xpvs=1./xinc
    c1xpvs=1.-xmin*c2xpvs
    do jx=1,nxpvs
      x=xmin+(jx-1)*xinc
      t=x
      tbpvs(jx)=fpvsx(t)
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This function computes saturation vapor pressure from the temperature.
!! A linear interpolation is done between values in a lookup table
!! computed in gpvs(). See documentation for fpvsx() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]  t     real, temperature in Kelvin
!\param[out] fpvs  real, saturation vapor pressure in Pascals
  elemental function fpvs(t)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvs         Compute saturation vapor pressure
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute saturation vapor pressure from the temperature.
!   A linear interpolation is done between values in a lookup table
!   computed in gpvs. See documentation for fpvsx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is almost 6 decimal places.
!   On the Cray, fpvs is about 4 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   pvs=fpvs(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvs       Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpvs
    real(krealfp),intent(in):: t
    integer jx
    real(krealfp) xj
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xpvs+c2xpvs*t,1._krealfp),real(nxpvs,krealfp))
    jx=min(xj,nxpvs-1._krealfp)
    fpvs=tbpvs(jx)+(xj-jx)*(tbpvs(jx+1)-tbpvs(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function computes saturation vapor pressure from the temperature.
!! A quadratic interpolation is done between values in a lookup table
!! computed in gpvs(). See documentation for fpvsx() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]  t        real, temperatue in Kelvin
!\param[out] fpvsq    real, saturation vapor pressure in Pascals
  elemental function fpvsq(t)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvsq        Compute saturation vapor pressure
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute saturation vapor pressure from the temperature.
!   A quadratic interpolation is done between values in a lookup table
!   computed in gpvs. See documentation for fpvsx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is almost 9 decimal places.
!   On the Cray, fpvsq is about 3 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             quadratic interpolation
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   pvs=fpvsq(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvsq      Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpvsq
    real(krealfp),intent(in):: t
    integer jx
    real(krealfp) xj,dxj,fj1,fj2,fj3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xpvs+c2xpvs*t,1._krealfp),real(nxpvs,krealfp))
    jx=min(max(nint(xj),2),nxpvs-1)
    dxj=xj-jx
    fj1=tbpvs(jx-1)
    fj2=tbpvs(jx)
    fj3=tbpvs(jx+1)
    fpvsq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function exactly computes saturation vapor pressure from temperature.
!! The saturation vapor pressure over either liquid and ice is computed
!! over liquid for temperatures above the triple point,
!! over ice for temperatures 20 degress below the triple point,
!! and a linear combination of the two for temperatures in between.
!! The water model assumes a perfect gas, constant specific heats
!! for gas, liquid and ice, and neglects the volume of the condensate.
!! The model does account for the variation of the latent heat
!! of condensation and sublimation with temperature.
!! The Clausius-Clapeyron equation is integrated from the triple point
!! to get the formula:
!!\n     pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
!!\n where tr is ttp/t and other values are physical constants.
!! The reference for this computation is Emanuel(1994), pages 116-117.
!! This function should be expanded inline in the calling routine.
!!\param[in]  t        real, temperature in Kelvin
!\param[out] fpvsx    real, saturation vapor pressure in Pascals
  elemental function fpvsx(t)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpvsx        Compute saturation vapor pressure
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Exactly compute saturation vapor pressure from temperature.
!   The saturation vapor pressure over either liquid and ice is computed
!   over liquid for temperatures above the triple point,
!   over ice for temperatures 20 degress below the triple point,
!   and a linear combination of the two for temperatures in between.
!   The water model assumes a perfect gas, constant specific heats
!   for gas, liquid and ice, and neglects the volume of the condensate.
!   The model does account for the variation of the latent heat
!   of condensation and sublimation with temperature.
!   The Clausius-Clapeyron equation is integrated from the triple point
!   to get the formula
!       pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
!   where tr is ttp/t and other values are physical constants.
!   The reference for this computation is Emanuel(1994), pages 116-117.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   pvs=fpvsx(t)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!
!   Output argument list:
!     fpvsx      Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpvsx
    real(krealfp),intent(in):: t
    real(krealfp),parameter:: tliq=con_ttp
    real(krealfp),parameter:: tice=con_ttp-20.0
    real(krealfp),parameter:: dldtl=con_cvap-con_cliq
    real(krealfp),parameter:: heatl=con_hvap
    real(krealfp),parameter:: xponal=-dldtl/con_rv
    real(krealfp),parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp)
    real(krealfp),parameter:: dldti=con_cvap-con_csol
    real(krealfp),parameter:: heati=con_hvap+con_hfus
    real(krealfp),parameter:: xponai=-dldti/con_rv
    real(krealfp),parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp)
    real(krealfp) tr,w,pvl,pvi
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    tr=con_ttp/t
    if(t.ge.tliq) then
      fpvsx=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
    elseif(t.lt.tice) then
      fpvsx=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
    else
      w=(t-tice)/(tliq-tice)
      pvl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
      pvi=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
      fpvsx=w*pvl+(1.-w)*pvi
    endif
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This subroutine computes dewpoint temperature table as a function of
!! vapor pressure for inlinable function ftdpl().
!! Exact dewpoint temperatures are calculated in subprogram ftdplxg().
!! The current implementation computes a table with a length
!! of 5001 for vapor pressures ranging from 1 to 10001 Pascals
!! giving a dewpoint temperature range of 208 to 319 Kelvin.
  subroutine gtdpl
!$$$     Subprogram Documentation Block
!
! Subprogram: gtdpl        Compute dewpoint temperature over liquid table
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute dewpoint temperature table as a function of
!   vapor pressure for inlinable function ftdpl.
!   Exact dewpoint temperatures are calculated in subprogram ftdplxg.
!   The current implementation computes a table with a length
!   of 5001 for vapor pressures ranging from 1 to 10001 Pascals
!   giving a dewpoint temperature range of 208 to 319 Kelvin.
!
! Program History Log:
!   91-05-07  Iredell
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
!
! Usage:  call gtdpl
!
! Subprograms called:
!   (ftdplxg)  inlinable function to compute dewpoint temperature over liquid
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx
    real(krealfp) xmin,xmax,xinc,t,x,pv
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=1
    xmax=10001
    xinc=(xmax-xmin)/(nxtdpl-1)
    c1xtdpl=1.-xmin/xinc
    c2xtdpl=1./xinc
    t=208.0
    do jx=1,nxtdpl
      x=xmin+(jx-1)*xinc
      pv=x
      t=ftdplxg(t,pv)
      tbtdpl(jx)=t
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This function Compute dewpoint temperature from vapor pressure.
!! A linear interpolation is done between values in a lookup table
!! computed in gtdpl(). See documentation for ftdplxg() for details.
!! Input values outside table range are reset to table extrema.
  elemental function ftdpl(pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdpl        Compute dewpoint temperature over liquid
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute dewpoint temperature from vapor pressure.
!   A linear interpolation is done between values in a lookup table
!   computed in gtdpl. See documentation for ftdplxg for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is better than 0.0005 Kelvin
!   for dewpoint temperatures greater than 250 Kelvin,
!   but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin.
!   On the Cray, ftdpl is about 75 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
!
! Usage:   tdpl=ftdpl(pv)
!
!   Input argument list:
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdpl      Real(krealfp) dewpoint temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdpl
    real(krealfp),intent(in):: pv
    integer jx
    real(krealfp) xj
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xtdpl+c2xtdpl*pv,1._krealfp),real(nxtdpl,krealfp))
    jx=min(xj,nxtdpl-1._krealfp)
    ftdpl=tbtdpl(jx)+(xj-jx)*(tbtdpl(jx+1)-tbtdpl(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function computes dewpoint temperature from vapor pressure.
!! A quadratic interpolation is done between values in a lookup table
!! computed in gtdpl(). See documentation for ftdplxg() for details.
!! Input values outside table range are reset to table extrema.
!!\param[in]  pv      real, vapor pressure in Pascals
!\param[out] ftdplq  real, dewpoint temperature in Kelvin
  elemental function ftdplq(pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdplq       Compute dewpoint temperature over liquid
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute dewpoint temperature from vapor pressure.
!   A quadratic interpolation is done between values in a lookup table
!   computed in gtdpl. see documentation for ftdplxg for details.
!   Input values outside table range are reset to table extrema.
!   the interpolation accuracy is better than 0.00001 Kelvin
!   for dewpoint temperatures greater than 250 Kelvin,
!   but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin.
!   On the Cray, ftdplq is about 60 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             quadratic interpolation
! 1999-03-01  Iredell             f90 module
!
! Usage:   tdpl=ftdplq(pv)
!
!   Input argument list:
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdplq     Real(krealfp) dewpoint temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdplq
    real(krealfp),intent(in):: pv
    integer jx
    real(krealfp) xj,dxj,fj1,fj2,fj3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xtdpl+c2xtdpl*pv,1._krealfp),real(nxtdpl,krealfp))
    jx=min(max(nint(xj),2),nxtdpl-1)
    dxj=xj-jx
    fj1=tbtdpl(jx-1)
    fj2=tbtdpl(jx)
    fj3=tbtdpl(jx+1)
    ftdplq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function exactly compute dewpoint temperature from vapor pressure.
!! An approximate dewpoint temperature for function ftdplxg()
!! is obtained using ftdpl() so gtdpl() must be already called.
!! See documentation for ftdplxg() for details.
!!\param[in]  pv      real, vapor pressure in Pascals
!\param[out] ftdplx  real, dewpoint temperature in Kelvin
  elemental function ftdplx(pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdplx       Compute dewpoint temperature over liquid
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: exactly compute dewpoint temperature from vapor pressure.
!   An approximate dewpoint temperature for function ftdplxg
!   is obtained using ftdpl so gtdpl must be already called.
!   See documentation for ftdplxg for details.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
!
! Usage:   tdpl=ftdplx(pv)
!
!   Input argument list:
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdplx     Real(krealfp) dewpoint temperature in Kelvin
!
! Subprograms called:
!   (ftdpl)    inlinable function to compute dewpoint temperature over liquid
!   (ftdplxg)  inlinable function to compute dewpoint temperature over liquid
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdplx
    real(krealfp),intent(in):: pv
    real(krealfp) tg
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    tg=ftdpl(pv)
    ftdplx=ftdplxg(tg,pv)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function exactly computes dewpoint temperature from vapor pressure.
!! A guess dewpoint temperature must be provided.
!! The water model assumes a perfect gas, constant specific heats
!! for gas and liquid, and neglects the volume of the liquid.
!! The model does account for the variation of the latent heat
!! of condensation with temperature.  The ice option is not included.
!! The Clausius-Clapeyron equation is integrated from the triple point
!! to get the formula:
!!\n  pvs=con_psat*(tr**xa)*exp(xb*(1.-tr))
!!\n  where tr is ttp/t and other values are physical constants.
!!\n The formula is inverted by iterating Newtonian approximations
!! for each pvs until t is found to within 1.e-6 Kelvin.
!! This function can be expanded inline in the calling routine.
!!\param[in]  tg       real, guess dewpoint temperature in Kelvin
!!\param[in]  pv       real, vapor pressure in Pascals
!\param[out] ftdplxg  real, dewpoint temperature in Kelvin
  elemental function ftdplxg(tg,pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdplxg      Compute dewpoint temperature over liquid
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Exactly compute dewpoint temperature from vapor pressure.
!   A guess dewpoint temperature must be provided.
!   The water model assumes a perfect gas, constant specific heats
!   for gas and liquid, and neglects the volume of the liquid.
!   The model does account for the variation of the latent heat
!   of condensation with temperature.  The ice option is not included.
!   The Clausius-Clapeyron equation is integrated from the triple point
!   to get the formula
!       pvs=con_psat*(tr**xa)*exp(xb*(1.-tr))
!   where tr is ttp/t and other values are physical constants.
!   The formula is inverted by iterating Newtonian approximations
!   for each pvs until t is found to within 1.e-6 Kelvin.
!   This function can be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
!
! Usage:   tdpl=ftdplxg(tg,pv)
!
!   Input argument list:
!     tg         Real(krealfp) guess dewpoint temperature in Kelvin
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdplxg    Real(krealfp) dewpoint temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdplxg
    real(krealfp),intent(in):: tg,pv
    real(krealfp),parameter:: terrm=1.e-6
    real(krealfp),parameter:: dldt=con_cvap-con_cliq
    real(krealfp),parameter:: heat=con_hvap
    real(krealfp),parameter:: xpona=-dldt/con_rv
    real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
    real(krealfp) t,tr,pvt,el,dpvt,terr
    integer i
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    t=tg
    do i=1,100
      tr=con_ttp/t
      pvt=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
      el=heat+dldt*(t-con_ttp)
      dpvt=el*pvt/(con_rv*t**2)
      terr=(pvt-pv)/dpvt
      t=t-terr
      if(abs(terr).le.terrm) exit
    enddo
    ftdplxg=t
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This subroutine computes dewpoint temperature table as a function of
!! vapor pressure for inlinable function ftdpi().
!! Exact dewpoint temperatures are calculated in subprogram ftdpixg().
!! The current implementation computes a table with a length
!! of 5001 for vapor pressures ranging from 0.1 to 1000.1 Pascals
!! giving a dewpoint temperature range of 197 to 279 Kelvin.
  subroutine gtdpi
!$$$     Subprogram Documentation Block
!
! Subprogram: gtdpi        Compute dewpoint temperature over ice table
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute dewpoint temperature table as a function of
!   vapor pressure for inlinable function ftdpi.
!   Exact dewpoint temperatures are calculated in subprogram ftdpixg.
!   The current implementation computes a table with a length
!   of 5001 for vapor pressures ranging from 0.1 to 1000.1 Pascals
!   giving a dewpoint temperature range of 197 to 279 Kelvin.
!
! Program History Log:
!   91-05-07  Iredell
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:  call gtdpi
!
! Subprograms called:
!   (ftdpixg)  inlinable function to compute dewpoint temperature over ice
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx
    real(krealfp) xmin,xmax,xinc,t,x,pv
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=0.1
    xmax=1000.1
    xinc=(xmax-xmin)/(nxtdpi-1)
    c1xtdpi=1.-xmin/xinc
    c2xtdpi=1./xinc
    t=197.0
    do jx=1,nxtdpi
      x=xmin+(jx-1)*xinc
      pv=x
      t=ftdpixg(t,pv)
      tbtdpi(jx)=t
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This subroutine computes dewpoint temperature from vapor pressure.
!! A linear interpolation is done between values in a lookup table
!! computed in gtdpi(). See documentation for ftdpixg for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]  pv     real, vapor pressure in Pascals
!\param[out] ftdpi  real, dewpoint temperature in Kelvin
  elemental function ftdpi(pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdpi        Compute dewpoint temperature over ice
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute dewpoint temperature from vapor pressure.
!   A linear interpolation is done between values in a lookup table
!   computed in gtdpi. See documentation for ftdpixg for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is better than 0.0005 Kelvin
!   for dewpoint temperatures greater than 250 Kelvin,
!   but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin.
!   On the Cray, ftdpi is about 75 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   tdpi=ftdpi(pv)
!
!   Input argument list:
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdpi      Real(krealfp) dewpoint temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdpi
    real(krealfp),intent(in):: pv
    integer jx
    real(krealfp) xj
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xtdpi+c2xtdpi*pv,1._krealfp),real(nxtdpi,krealfp))
    jx=min(xj,nxtdpi-1._krealfp)
    ftdpi=tbtdpi(jx)+(xj-jx)*(tbtdpi(jx+1)-tbtdpi(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function computes dewpoint temperature from vapor pressure.
!! A quadratic interpolation is done between values in a lookup table
!! computed in gtdpi(). see documentation for ftdpixg() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]  pv       real, vapor pressure in Pascals
!\param[out] ftdpiq   real, dewpoint temperature in Kelvin
  elemental function ftdpiq(pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdpiq       Compute dewpoint temperature over ice
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute dewpoint temperature from vapor pressure.
!   A quadratic interpolation is done between values in a lookup table
!   computed in gtdpi. see documentation for ftdpixg for details.
!   Input values outside table range are reset to table extrema.
!   the interpolation accuracy is better than 0.00001 Kelvin
!   for dewpoint temperatures greater than 250 Kelvin,
!   but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin.
!   On the Cray, ftdpiq is about 60 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             quadratic interpolation
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   tdpi=ftdpiq(pv)
!
!   Input argument list:
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdpiq     Real(krealfp) dewpoint temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdpiq
    real(krealfp),intent(in):: pv
    integer jx
    real(krealfp) xj,dxj,fj1,fj2,fj3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xtdpi+c2xtdpi*pv,1._krealfp),real(nxtdpi,krealfp))
    jx=min(max(nint(xj),2),nxtdpi-1)
    dxj=xj-jx
    fj1=tbtdpi(jx-1)
    fj2=tbtdpi(jx)
    fj3=tbtdpi(jx+1)
    ftdpiq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function exactly computes dewpoint temperature from vapor pressure.
!! An approximate dewpoint temperature for function ftdpixg()
!! is obtained using ftdpi() so gtdpi() must be already called.
!! See documentation for ftdpixg() for details.
!>\param[in]  pv   real, vapor pressure in Pascals
!\param[out] ftdpix   real, dewpoint temperature in Kelvin
  elemental function ftdpix(pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdpix       Compute dewpoint temperature over ice
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: exactly compute dewpoint temperature from vapor pressure.
!   An approximate dewpoint temperature for function ftdpixg
!   is obtained using ftdpi so gtdpi must be already called.
!   See documentation for ftdpixg for details.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   tdpi=ftdpix(pv)
!
!   Input argument list:
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdpix     Real(krealfp) dewpoint temperature in Kelvin
!
! Subprograms called:
!   (ftdpi)    inlinable function to compute dewpoint temperature over ice
!   (ftdpixg)  inlinable function to compute dewpoint temperature over ice
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdpix
    real(krealfp),intent(in):: pv
    real(krealfp) tg
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    tg=ftdpi(pv)
    ftdpix=ftdpixg(tg,pv)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function exactly computes dewpoint temperature from vapor pressure.
!! A guess dewpoint temperature must be provided.
!! The water model assumes a perfect gas, constant specific heats
!! for gas and ice, and neglects the volume of the ice.
!! The model does account for the variation of the latent heat
!! of sublimation with temperature.  The liquid option is not included.
!! The Clausius-Clapeyron equation is integrated from the triple point
!! to get the formula:
!!\n   pvs=con_psat*(tr**xa)*exp(xb*(1.-tr))
!!\n  where tr is ttp/t and other values are physical constants.
!! The formula is inverted by iterating Newtonian approximations
!! for each pvs until t is found to within 1.e-6 Kelvin.
!!   This function can be expanded inline in the calling routine.
!>\param[in]  tg   real, guess dewpoint temperature in Kelvin
!>\param[in]  pv   real, vapor pressure in Pascals
!\param[out] ftdpixg   real, dewpoint temperature in Kelvin
  elemental function ftdpixg(tg,pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdpixg      Compute dewpoint temperature over ice
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Exactly compute dewpoint temperature from vapor pressure.
!   A guess dewpoint temperature must be provided.
!   The water model assumes a perfect gas, constant specific heats
!   for gas and ice, and neglects the volume of the ice.
!   The model does account for the variation of the latent heat
!   of sublimation with temperature.  The liquid option is not included.
!   The Clausius-Clapeyron equation is integrated from the triple point
!   to get the formula
!       pvs=con_psat*(tr**xa)*exp(xb*(1.-tr))
!   where tr is ttp/t and other values are physical constants.
!   The formula is inverted by iterating Newtonian approximations
!   for each pvs until t is found to within 1.e-6 Kelvin.
!   This function can be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   tdpi=ftdpixg(tg,pv)
!
!   Input argument list:
!     tg         Real(krealfp) guess dewpoint temperature in Kelvin
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdpixg    Real(krealfp) dewpoint temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdpixg
    real(krealfp),intent(in):: tg,pv
    real(krealfp),parameter:: terrm=1.e-6
    real(krealfp),parameter:: dldt=con_cvap-con_csol
    real(krealfp),parameter:: heat=con_hvap+con_hfus
    real(krealfp),parameter:: xpona=-dldt/con_rv
    real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp)
    real(krealfp) t,tr,pvt,el,dpvt,terr
    integer i
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    t=tg
    do i=1,100
      tr=con_ttp/t
      pvt=con_psat*(tr**xpona)*exp(xponb*(1.-tr))
      el=heat+dldt*(t-con_ttp)
      dpvt=el*pvt/(con_rv*t**2)
      terr=(pvt-pv)/dpvt
      t=t-terr
      if(abs(terr).le.terrm) exit
    enddo
    ftdpixg=t
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This subroutine computes dewpoint temperature table as a function of
!! vapor pressure for inlinable function ftdp().
!! Exact dewpoint temperatures are calculated in subprogram ftdpxg().
!! The current implementation computes a table with a length
!! of 5001 for vapor pressures ranging from 0.5 to 1000.5 Pascals
!! giving a dewpoint temperature range of 208 to 319 Kelvin.
  subroutine gtdp
!$$$     Subprogram Documentation Block
!
! Subprogram: gtdp         Compute dewpoint temperature table
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute dewpoint temperature table as a function of
!   vapor pressure for inlinable function ftdp.
!   Exact dewpoint temperatures are calculated in subprogram ftdpxg.
!   The current implementation computes a table with a length
!   of 5001 for vapor pressures ranging from 0.5 to 1000.5 Pascals
!   giving a dewpoint temperature range of 208 to 319 Kelvin.
!
! Program History Log:
!   91-05-07  Iredell
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:  call gtdp
!
! Subprograms called:
!   (ftdpxg)   inlinable function to compute dewpoint temperature
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx
    real(krealfp) xmin,xmax,xinc,t,x,pv
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=0.5
    xmax=10000.5
    xinc=(xmax-xmin)/(nxtdp-1)
    c1xtdp=1.-xmin/xinc
    c2xtdp=1./xinc
    t=208.0
    do jx=1,nxtdp
      x=xmin+(jx-1)*xinc
      pv=x
      t=ftdpxg(t,pv)
      tbtdp(jx)=t
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This function computes dewpoint temperature from vapor pressure.
!! A linear interpolation is done between values in a lookup table
!! computed in gtdp(). See documentation for ftdpxg() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]  pv  real, vapor pressure in Pascals
!\param[out] ftdp  real, dewpoint temperature in Kelvin
  elemental function ftdp(pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdp         Compute dewpoint temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute dewpoint temperature from vapor pressure.
!   A linear interpolation is done between values in a lookup table
!   computed in gtdp. See documentation for ftdpxg for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is better than 0.0005 Kelvin
!   for dewpoint temperatures greater than 250 Kelvin,
!   but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin.
!   On the Cray, ftdp is about 75 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   tdp=ftdp(pv)
!
!   Input argument list:
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdp       Real(krealfp) dewpoint temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdp
    real(krealfp),intent(in):: pv
    integer jx
    real(krealfp) xj
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xtdp+c2xtdp*pv,1._krealfp),real(nxtdp,krealfp))
    jx=min(xj,nxtdp-1._krealfp)
    ftdp=tbtdp(jx)+(xj-jx)*(tbtdp(jx+1)-tbtdp(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function computes dewpoint temperature from vapor pressure.
!! A quadratic interpolation is done between values in a lookup table
!! computed in gtdp(). See documentation for ftdpxg() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]  pv   real, vapor pressure in Pascals
!\param[out] ftdpq   real, dewpoint temperature in Kelvin
  elemental function ftdpq(pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdpq        Compute dewpoint temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute dewpoint temperature from vapor pressure.
!   A quadratic interpolation is done between values in a lookup table
!   computed in gtdp. see documentation for ftdpxg for details.
!   Input values outside table range are reset to table extrema.
!   the interpolation accuracy is better than 0.00001 Kelvin
!   for dewpoint temperatures greater than 250 Kelvin,
!   but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin.
!   On the Cray, ftdpq is about 60 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             quadratic interpolation
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   tdp=ftdpq(pv)
!
!   Input argument list:
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdpq      Real(krealfp) dewpoint temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdpq
    real(krealfp),intent(in):: pv
    integer jx
    real(krealfp) xj,dxj,fj1,fj2,fj3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xtdp+c2xtdp*pv,1._krealfp),real(nxtdp,krealfp))
    jx=min(max(nint(xj),2),nxtdp-1)
    dxj=xj-jx
    fj1=tbtdp(jx-1)
    fj2=tbtdp(jx)
    fj3=tbtdp(jx+1)
    ftdpq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function exactly computes dewpoint temperature from vapor pressure.
!! An approximate dewpoint temperature for function ftdpxg()
!! is obtained using ftdp() so gtdp() must be already called.
!! See documentation for ftdpxg() for details.
!>\param[in]  pv   real, vapor pressure in Pascals
!\param[out] ftdpx    real, dewpoint temperature in Kelvin
  elemental function ftdpx(pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdpx        Compute dewpoint temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: exactly compute dewpoint temperature from vapor pressure.
!   An approximate dewpoint temperature for function ftdpxg
!   is obtained using ftdp so gtdp must be already called.
!   See documentation for ftdpxg for details.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   tdp=ftdpx(pv)
!
!   Input argument list:
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdpx      Real(krealfp) dewpoint temperature in Kelvin
!
! Subprograms called:
!   (ftdp)     inlinable function to compute dewpoint temperature
!   (ftdpxg)   inlinable function to compute dewpoint temperature
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdpx
    real(krealfp),intent(in):: pv
    real(krealfp) tg
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    tg=ftdp(pv)
    ftdpx=ftdpxg(tg,pv)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function exactly computes dewpoint temperature from vapor pressure.
!! A guess dewpoint temperature must be provided.
!! The saturation vapor pressure over either liquid and ice is computed
!! over liquid for temperatures above the triple point,
!! over ice for temperatures 20 degress below the triple point,
!! and a linear combination of the two for temperatures in between.
!! The water model assumes a perfect gas, constant specific heats
!! for gas, liquid and ice, and neglects the volume of the condensate.
!! The model does account for the variation of the latent heat
!! of condensation and sublimation with temperature.
!! The Clausius-Clapeyron equation is integrated from the triple point
!! to get the formula:
!!\n   pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
!!\n where tr is ttp/t and other values are physical constants.
!!\n The reference for this decision is Emanuel(1994), pages 116-117.
!! The formula is inverted by iterating Newtonian approximations
!! for each pvs until t is found to within 1.e-6 Kelvin.
!! This function can be expanded inline in the calling routine.
!>\param[in]  tg   real, guess dewpoint temperature in Kelvin
!>\param[in]  pv   real, vapor pressure in Pascals
!\param[out] ftdpxg    real, dewpoint temperature in Kelvin
  elemental function ftdpxg(tg,pv)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftdpxg       Compute dewpoint temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Exactly compute dewpoint temperature from vapor pressure.
!   A guess dewpoint temperature must be provided.
!   The saturation vapor pressure over either liquid and ice is computed
!   over liquid for temperatures above the triple point,
!   over ice for temperatures 20 degress below the triple point,
!   and a linear combination of the two for temperatures in between.
!   The water model assumes a perfect gas, constant specific heats
!   for gas, liquid and ice, and neglects the volume of the condensate.
!   The model does account for the variation of the latent heat
!   of condensation and sublimation with temperature.
!   The Clausius-Clapeyron equation is integrated from the triple point
!   to get the formula
!       pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr))
!   where tr is ttp/t and other values are physical constants.
!   The reference for this decision is Emanuel(1994), pages 116-117.
!   The formula is inverted by iterating Newtonian approximations
!   for each pvs until t is found to within 1.e-6 Kelvin.
!   This function can be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
! 2001-02-26  Iredell             ice phase
!
! Usage:   tdp=ftdpxg(tg,pv)
!
!   Input argument list:
!     tg         Real(krealfp) guess dewpoint temperature in Kelvin
!     pv         Real(krealfp) vapor pressure in Pascals
!
!   Output argument list:
!     ftdpxg     Real(krealfp) dewpoint temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftdpxg
    real(krealfp),intent(in):: tg,pv
    real(krealfp),parameter:: terrm=1.e-6
    real(krealfp),parameter:: tliq=con_ttp
    real(krealfp),parameter:: tice=con_ttp-20.0
    real(krealfp),parameter:: dldtl=con_cvap-con_cliq
    real(krealfp),parameter:: heatl=con_hvap
    real(krealfp),parameter:: xponal=-dldtl/con_rv
    real(krealfp),parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp)
    real(krealfp),parameter:: dldti=con_cvap-con_csol
    real(krealfp),parameter:: heati=con_hvap+con_hfus
    real(krealfp),parameter:: xponai=-dldti/con_rv
    real(krealfp),parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp)
    real(krealfp) t,tr,w,pvtl,pvti,pvt,ell,eli,el,dpvt,terr
    integer i
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    t=tg
    do i=1,100
      tr=con_ttp/t
      if(t.ge.tliq) then
        pvt=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
        el=heatl+dldtl*(t-con_ttp)
        dpvt=el*pvt/(con_rv*t**2)
      elseif(t.lt.tice) then
        pvt=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
        el=heati+dldti*(t-con_ttp)
        dpvt=el*pvt/(con_rv*t**2)
      else
        w=(t-tice)/(tliq-tice)
        pvtl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
        pvti=con_psat*(tr**xponai)*exp(xponbi*(1.-tr))
        pvt=w*pvtl+(1.-w)*pvti
        ell=heatl+dldtl*(t-con_ttp)
        eli=heati+dldti*(t-con_ttp)
        dpvt=(w*ell*pvtl+(1.-w)*eli*pvti)/(con_rv*t**2)
      endif
      terr=(pvt-pv)/dpvt
      t=t-terr
      if(abs(terr).le.terrm) exit
    enddo
    ftdpxg=t
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This subroutine computes equivalent potential temperature table
!! as a function of LCL temperature and pressure over 1e5 Pa
!! to the kappa power for function fthe().
!! Equivalent potential temperatures are calculated in subprogram fthex()
!! the current implementation computes a table with a first dimension
!! of 241 for temperatures ranging from 183.16 to 303.16 Kelvin
!! and a second dimension of 151 for pressure over 1e5 Pa
!! to the kappa power ranging from \f$0.04^{rocp}\f$ to \f$1.10^{rocp}\f$.
  subroutine gthe
!$$$     Subprogram Documentation Block
!
! Subprogram: gthe        Compute equivalent potential temperature table
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute equivalent potential temperature table
!   as a function of LCL temperature and pressure over 1e5 Pa
!   to the kappa power for function fthe.
!   Equivalent potential temperatures are calculated in subprogram fthex
!   the current implementation computes a table with a first dimension
!   of 241 for temperatures ranging from 183.16 to 303.16 Kelvin
!   and a second dimension of 151 for pressure over 1e5 Pa
!   to the kappa power ranging from 0.04**rocp to 1.10**rocp.
!
! Program History Log:
!   91-05-07  Iredell
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
!
! Usage:  call gthe
!
! Subprograms called:
!   (fthex)    inlinable function to compute equiv. pot. temperature
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx,jy
    real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,pk,t
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=con_ttp-90._krealfp
    xmax=con_ttp+30._krealfp
    ymin=0.04_krealfp**con_rocp
    ymax=1.10_krealfp**con_rocp
    xinc=(xmax-xmin)/(nxthe-1)
    c1xthe=1.-xmin/xinc
    c2xthe=1./xinc
    yinc=(ymax-ymin)/(nythe-1)
    c1ythe=1.-ymin/yinc
    c2ythe=1./yinc
    do jy=1,nythe
      y=ymin+(jy-1)*yinc
      pk=y
      do jx=1,nxthe
        x=xmin+(jx-1)*xinc
        t=x
        tbthe(jx,jy)=fthex(t,pk)
      enddo
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This function computes equivalent potential temperature at the LCL
!! from temperature and pressure over 1e5 Pa to the kappa power.
!! A bilinear interpolation is done between values in a lookup table
!! computed in gthe(). see documentation for fthex() for details.
!! Input values outside table range are reset to table extrema,
!!   except zero is returned for too cold or high LCLs.
!>\param[in]  t  real, LCL temperature in Kelvin
!>\param[in]  pk  real, LCL pressure over 1e5 Pa to the kappa power
!\param[out] fthe     real, equivalent potential temperature in Kelvin
  elemental function fthe(t,pk)
!$$$     Subprogram Documentation Block
!
! Subprogram: fthe         Compute equivalent potential temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute equivalent potential temperature at the LCL
!   from temperature and pressure over 1e5 Pa to the kappa power.
!   A bilinear interpolation is done between values in a lookup table
!   computed in gthe. see documentation for fthex for details.
!   Input values outside table range are reset to table extrema,
!   except zero is returned for too cold or high LCLs.
!   The interpolation accuracy is better than 0.01 Kelvin.
!   On the Cray, fthe is almost 6 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
!
! Usage:   the=fthe(t,pk)
!
!   Input argument list:
!     t          Real(krealfp) LCL temperature in Kelvin
!     pk         Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
!
!   Output argument list:
!     fthe       Real(krealfp) equivalent potential temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fthe
    real(krealfp),intent(in):: t,pk
    integer jx,jy
    real(krealfp) xj,yj,ftx1,ftx2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(c1xthe+c2xthe*t,real(nxthe,krealfp))
    yj=min(c1ythe+c2ythe*pk,real(nythe,krealfp))
    if(xj.ge.1..and.yj.ge.1.) then
      jx=min(xj,nxthe-1._krealfp)
      jy=min(yj,nythe-1._krealfp)
      ftx1=tbthe(jx,jy)+(xj-jx)*(tbthe(jx+1,jy)-tbthe(jx,jy))
      ftx2=tbthe(jx,jy+1)+(xj-jx)*(tbthe(jx+1,jy+1)-tbthe(jx,jy+1))
      fthe=ftx1+(yj-jy)*(ftx2-ftx1)
    else
      fthe=0.
    endif
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function computes equivalent potential temperature at the LCL
!! from temperature and pressure over 1e5 Pa to the kappa power.
!! A biquadratic interpolation is done between values in a lookup table
!! computed in gthe(). see documentation for fthex() for details.
!! Input values outside table range are reset to table extrema,
!! except zero is returned for too cold or high LCLs.
!>\param[in]  t    real, LCL temperature in Kelvin
!>\param[in]  pk   real, LCL pressure over 1e5 Pa to the kappa power
!\param[out] ftheq       real, equivalent potential temperature in Kelvin
  elemental function ftheq(t,pk)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftheq        Compute equivalent potential temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute equivalent potential temperature at the LCL
!   from temperature and pressure over 1e5 Pa to the kappa power.
!   A biquadratic interpolation is done between values in a lookup table
!   computed in gthe. see documentation for fthex for details.
!   Input values outside table range are reset to table extrema,
!   except zero is returned for too cold or high LCLs.
!   The interpolation accuracy is better than 0.0002 Kelvin.
!   On the Cray, ftheq is almost 3 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             quadratic interpolation
! 1999-03-01  Iredell             f90 module
!
! Usage:   the=ftheq(t,pk)
!
!   Input argument list:
!     t          Real(krealfp) LCL temperature in Kelvin
!     pk         Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
!
!   Output argument list:
!     ftheq      Real(krealfp) equivalent potential temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftheq
    real(krealfp),intent(in):: t,pk
    integer jx,jy
    real(krealfp) xj,yj,dxj,dyj
    real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33
    real(krealfp) ftx1,ftx2,ftx3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(c1xthe+c2xthe*t,real(nxthe,krealfp))
    yj=min(c1ythe+c2ythe*pk,real(nythe,krealfp))
    if(xj.ge.1..and.yj.ge.1.) then
      jx=min(max(nint(xj),2),nxthe-1)
      jy=min(max(nint(yj),2),nythe-1)
      dxj=xj-jx
      dyj=yj-jy
      ft11=tbthe(jx-1,jy-1)
      ft12=tbthe(jx-1,jy)
      ft13=tbthe(jx-1,jy+1)
      ft21=tbthe(jx,jy-1)
      ft22=tbthe(jx,jy)
      ft23=tbthe(jx,jy+1)
      ft31=tbthe(jx+1,jy-1)
      ft32=tbthe(jx+1,jy)
      ft33=tbthe(jx+1,jy+1)
      ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21
      ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22
      ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23
      ftheq=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2
    else
      ftheq=0.
    endif
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
! elemental function fthex(t,pk)
!> This function exactly computes equivalent potential temperature at the LCL
!! from temperature and pressure over 1e5 Pa to the kappa power.
!! Equivalent potential temperature is constant for a saturated parcel
!! rising adiabatically up a moist adiabat when the heat and mass
!! of the condensed water are neglected.  Ice is also neglected.
!! The formula for equivalent potential temperature (Holton) is
!!\n  the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd))
!!\n where t is the temperature, pv is the saturated vapor pressure,
!! pd is the dry pressure p-pv, el is the temperature dependent
!! latent heat of condensation hvap+dldt*(t-ttp), and other values
!! are physical constants defined in parameter statements in the code.
!! Zero is returned if the input values make saturation impossible.
!! This function should be expanded inline in the calling routine.
!>\param[in]  t  real, LCL temperature in Kelvin
!>\param[in]  pk  real, LCL pressure over 1e5 Pa to the kappa power
!\param[out] fthex     real, equivalent potential temperature in Kelvin
            function fthex(t,pk)
!$$$     Subprogram Documentation Block
!
! Subprogram: fthex        Compute equivalent potential temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Exactly compute equivalent potential temperature at the LCL
!   from temperature and pressure over 1e5 Pa to the kappa power.
!   Equivalent potential temperature is constant for a saturated parcel
!   rising adiabatically up a moist adiabat when the heat and mass
!   of the condensed water are neglected.  Ice is also neglected.
!   The formula for equivalent potential temperature (Holton) is
!       the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd))
!   where t is the temperature, pv is the saturated vapor pressure,
!   pd is the dry pressure p-pv, el is the temperature dependent
!   latent heat of condensation hvap+dldt*(t-ttp), and other values
!   are physical constants defined in parameter statements in the code.
!   Zero is returned if the input values make saturation impossible.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
!
! Usage:   the=fthex(t,pk)
!
!   Input argument list:
!     t          Real(krealfp) LCL temperature in Kelvin
!     pk         Real(krealfp) LCL pressure over 1e5 Pa to the kappa power
!
!   Output argument list:
!     fthex      Real(krealfp) equivalent potential temperature in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fthex
    real(krealfp),intent(in):: t,pk
    real(krealfp) p,tr,pv,pd,el,expo,expmax
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    p=pk**con_cpor
    tr=con_ttp/t
    pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr))
    pd=p-pv
    if(pd.gt.pv) then
      el=con_hvap+con_dldt*(t-con_ttp)
      expo=el*con_eps*pv/(con_cp*t*pd)
      fthex=t*pd**(-con_rocp)*exp(expo)
    else
      fthex=0.
    endif
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This subroutine computes temperature and specific humidity tables
!! as a function of equivalent potential temperature and
!! pressure over 1e5 Pa to the kappa power for subprogram stma().
!! Exact parcel temperatures are calculated in subprogram stmaxg().
!! The current implementation computes a table with a first dimension
!! of 151 for equivalent potential temperatures ranging from 200 to 500
!! Kelvin and a second dimension of 121 for pressure over 1e5 Pa
!! to the kappa power ranging from \f$0.01^{rocp}\f$ to \f$1.10^{rocp}\f$.
  subroutine gtma
!$$$     Subprogram Documentation Block
!
! Subprogram: gtma         Compute moist adiabat tables
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute temperature and specific humidity tables
!   as a function of equivalent potential temperature and
!   pressure over 1e5 Pa to the kappa power for subprogram stma.
!   Exact parcel temperatures are calculated in subprogram stmaxg.
!   The current implementation computes a table with a first dimension
!   of 151 for equivalent potential temperatures ranging from 200 to 500
!   Kelvin and a second dimension of 121 for pressure over 1e5 Pa
!   to the kappa power ranging from 0.01**rocp to 1.10**rocp.
!
! Program History Log:
!   91-05-07  Iredell
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
!
! Usage:  call gtma
!
! Subprograms called:
!   (stmaxg)   inlinable subprogram to compute parcel temperature
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx,jy
    real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,pk,the,t,q,tg
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=200._krealfp
    xmax=500._krealfp
    ymin=0.01_krealfp**con_rocp
    ymax=1.10_krealfp**con_rocp
    xinc=(xmax-xmin)/(nxma-1)
    c1xma=1.-xmin/xinc
    c2xma=1./xinc
    yinc=(ymax-ymin)/(nyma-1)
    c1yma=1.-ymin/yinc
    c2yma=1./yinc
    do jy=1,nyma
      y=ymin+(jy-1)*yinc
      pk=y
      tg=xmin*y
      do jx=1,nxma
        x=xmin+(jx-1)*xinc
        the=x
        call stmaxg(tg,the,pk,t,q)
        tbtma(jx,jy)=t
        tbqma(jx,jy)=q
        tg=t
      enddo
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This subroutine computes temperature and specific humidity of a parcel
!! lifted up a moist adiabat from equivalent potential temperature
!! at the LCL and pressure over 1e5 Pa to the kappa power.
!! Bilinear interpolations are done between values in a lookup table
!! computed in gtma(). See documentation for stmaxg() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]  the    real, equivalent potential temperature in Kelvin
!>\param[in]  pk     real, pressure over 1e5 Pa to the kappa power
!>\param[out] tma    real, parcel temperature in Kelvin
!>\param[out] qma    real, parcel specific humidity in kg/kg
  elemental subroutine stma(the,pk,tma,qma)
!$$$     Subprogram Documentation Block
!
! Subprogram: stma         Compute moist adiabat temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute temperature and specific humidity of a parcel
!   lifted up a moist adiabat from equivalent potential temperature
!   at the LCL and pressure over 1e5 Pa to the kappa power.
!   Bilinear interpolations are done between values in a lookup table
!   computed in gtma. See documentation for stmaxg for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is better than 0.01 Kelvin
!   and 5.e-6 kg/kg for temperature and humidity, respectively.
!   On the Cray, stma is about 35 times faster than exact calculation.
!   This subprogram should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             expand table
! 1999-03-01  Iredell             f90 module
!
! Usage:  call stma(the,pk,tma,qma)
!
!   Input argument list:
!     the        Real(krealfp) equivalent potential temperature in Kelvin
!     pk         Real(krealfp) pressure over 1e5 Pa to the kappa power
!
!   Output argument list:
!     tma        Real(krealfp) parcel temperature in Kelvin
!     qma        Real(krealfp) parcel specific humidity in kg/kg
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp),intent(in):: the,pk
    real(krealfp),intent(out):: tma,qma
    integer jx,jy
    real(krealfp) xj,yj,ftx1,ftx2,qx1,qx2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xma+c2xma*the,1._krealfp),real(nxma,krealfp))
    yj=min(max(c1yma+c2yma*pk,1._krealfp),real(nyma,krealfp))
    jx=min(xj,nxma-1._krealfp)
    jy=min(yj,nyma-1._krealfp)
    ftx1=tbtma(jx,jy)+(xj-jx)*(tbtma(jx+1,jy)-tbtma(jx,jy))
    ftx2=tbtma(jx,jy+1)+(xj-jx)*(tbtma(jx+1,jy+1)-tbtma(jx,jy+1))
    tma=ftx1+(yj-jy)*(ftx2-ftx1)
    qx1=tbqma(jx,jy)+(xj-jx)*(tbqma(jx+1,jy)-tbqma(jx,jy))
    qx2=tbqma(jx,jy+1)+(xj-jx)*(tbqma(jx+1,jy+1)-tbqma(jx,jy+1))
    qma=qx1+(yj-jy)*(qx2-qx1)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This subroutine computes temperature and specific humidity of a parcel
!! lifted up a moist adiabat from equivalent potential temperature
!! at the LCL and pressure over 1e5 Pa to the kappa power.
!! Biquadratic interpolations are done between values in a lookup table
!! computed in gtma(). See documentation for stmaxg() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]  the    real, equivalent potential temperature in Kelvin
!>\param[in]  pk     real, pressure over 1e5 Pa to the kappa power
!>\param[out] tma   real, parcel temperature in Kelvin
!>\param[out] qma    real, parcel specific humidity in kg/kg
  elemental subroutine stmaq(the,pk,tma,qma)
!$$$     Subprogram Documentation Block
!
! Subprogram: stmaq        Compute moist adiabat temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute temperature and specific humidity of a parcel
!   lifted up a moist adiabat from equivalent potential temperature
!   at the LCL and pressure over 1e5 Pa to the kappa power.
!   Biquadratic interpolations are done between values in a lookup table
!   computed in gtma. See documentation for stmaxg for details.
!   Input values outside table range are reset to table extrema.
!   the interpolation accuracy is better than 0.0005 Kelvin
!   and 1.e-7 kg/kg for temperature and humidity, respectively.
!   On the Cray, stmaq is about 25 times faster than exact calculation.
!   This subprogram should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             quadratic interpolation
! 1999-03-01  Iredell             f90 module
!
! Usage:  call stmaq(the,pk,tma,qma)
!
!   Input argument list:
!     the        Real(krealfp) equivalent potential temperature in Kelvin
!     pk         Real(krealfp) pressure over 1e5 Pa to the kappa power
!
!   Output argument list:
!     tmaq       Real(krealfp) parcel temperature in Kelvin
!     qma        Real(krealfp) parcel specific humidity in kg/kg
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp),intent(in):: the,pk
    real(krealfp),intent(out):: tma,qma
    integer jx,jy
    real(krealfp) xj,yj,dxj,dyj
    real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33
    real(krealfp) ftx1,ftx2,ftx3
    real(krealfp) q11,q12,q13,q21,q22,q23,q31,q32,q33,qx1,qx2,qx3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xma+c2xma*the,1._krealfp),real(nxma,krealfp))
    yj=min(max(c1yma+c2yma*pk,1._krealfp),real(nyma,krealfp))
    jx=min(max(nint(xj),2),nxma-1)
    jy=min(max(nint(yj),2),nyma-1)
    dxj=xj-jx
    dyj=yj-jy
    ft11=tbtma(jx-1,jy-1)
    ft12=tbtma(jx-1,jy)
    ft13=tbtma(jx-1,jy+1)
    ft21=tbtma(jx,jy-1)
    ft22=tbtma(jx,jy)
    ft23=tbtma(jx,jy+1)
    ft31=tbtma(jx+1,jy-1)
    ft32=tbtma(jx+1,jy)
    ft33=tbtma(jx+1,jy+1)
    ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21
    ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22
    ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23
    tma=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2
    q11=tbqma(jx-1,jy-1)
    q12=tbqma(jx-1,jy)
    q13=tbqma(jx-1,jy+1)
    q21=tbqma(jx,jy-1)
    q22=tbqma(jx,jy)
    q23=tbqma(jx,jy+1)
    q31=tbqma(jx+1,jy-1)
    q32=tbqma(jx+1,jy)
    q33=tbqma(jx+1,jy+1)
    qx1=(((q31+q11)/2-q21)*dxj+(q31-q11)/2)*dxj+q21
    qx2=(((q32+q12)/2-q22)*dxj+(q32-q12)/2)*dxj+q22
    qx3=(((q33+q13)/2-q23)*dxj+(q33-q13)/2)*dxj+q23
    qma=(((qx3+qx1)/2-qx2)*dyj+(qx3-qx1)/2)*dyj+qx2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This subroutine exactly computes temperature and humidity of a parcel
!! lifted up a moist adiabat from equivalent potential temperature
!! at the LCL and pressure over 1e5 Pa to the kappa power.
!! An approximate parcel temperature for subprogram stmaxg()
!! is obtained using stma() so gtma() must be already called.
!! See documentation for stmaxg() for details.
!>\param[in]  the   real, equivalent potential temperature in Kelvin
!>\param[in]  pk    real, pressure over 1e5 Pa to the kappa power
!>\param[out] tma   real, parcel temperature in Kelvin
!>\param[out] qma   real, parcel specific humidity in kg/kg
  subroutine stmax(the,pk,tma,qma)
!$$$     Subprogram Documentation Block
!
! Subprogram: stmax        Compute moist adiabat temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Exactly compute temperature and humidity of a parcel
!   lifted up a moist adiabat from equivalent potential temperature
!   at the LCL and pressure over 1e5 Pa to the kappa power.
!   An approximate parcel temperature for subprogram stmaxg
!   is obtained using stma so gtma must be already called.
!   See documentation for stmaxg for details.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
!
! Usage:  call stmax(the,pk,tma,qma)
!
!   Input argument list:
!     the        Real(krealfp) equivalent potential temperature in Kelvin
!     pk         Real(krealfp) pressure over 1e5 Pa to the kappa power
!
!   Output argument list:
!     tma        Real(krealfp) parcel temperature in Kelvin
!     qma        Real(krealfp) parcel specific humidity in kg/kg
!
! Subprograms called:
!   (stma)     inlinable subprogram to compute parcel temperature
!   (stmaxg)   inlinable subprogram to compute parcel temperature
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp),intent(in):: the,pk
    real(krealfp),intent(out):: tma,qma
    real(krealfp) tg,qg
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    call stma(the,pk,tg,qg)
    call stmaxg(tg,the,pk,tma,qma)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This subroutine exactly computes temperature and humidity of a parcel
!! lifted up a moist adiabat from equivalent potential temperature
!! at the LCL and pressure over 1e5 Pa to the kappa power.
!! A guess parcel temperature must be provided.
!! Equivalent potential temperature is constant for a saturated parcel
!! rising adiabatically up a moist adiabat when the heat and mass
!! of the condensed water are neglected.  Ice is also neglected.
!! The formula for equivalent potential temperature (Holton) is
!!\n  the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd))
!!\n where t is the temperature, pv is the saturated vapor pressure,
!! pd is the dry pressure p-pv, el is the temperature dependent
!! latent heat of condensation hvap+dldt*(t-ttp), and other values
!! are physical constants defined in parameter statements in the code.
!! The formula is inverted by iterating Newtonian approximations
!! for each the and p until t is found to within 1.e-4 Kelvin.
!! The specific humidity is then computed from pv and pd.
!! This subprogram can be expanded inline in the calling routine.
!>\param[in]  tg   real, guess parcel temperature in Kelvin
!>\param[in]  the  real, equivalent potential temperature in Kelvin
!>\param[in]  pk   real, pressure over 1e5 Pa to the kappa power
!>\param[out] tma  real, parcel temperature in Kelvin 
!>\param[out] qma  real, parcel specific humidity in kg/kg
  subroutine stmaxg(tg,the,pk,tma,qma)
!$$$     Subprogram Documentation Block
!
! Subprogram: stmaxg       Compute moist adiabat temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: exactly compute temperature and humidity of a parcel
!   lifted up a moist adiabat from equivalent potential temperature
!   at the LCL and pressure over 1e5 Pa to the kappa power.
!   A guess parcel temperature must be provided.
!   Equivalent potential temperature is constant for a saturated parcel
!   rising adiabatically up a moist adiabat when the heat and mass
!   of the condensed water are neglected.  Ice is also neglected.
!   The formula for equivalent potential temperature (Holton) is
!       the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd))
!   where t is the temperature, pv is the saturated vapor pressure,
!   pd is the dry pressure p-pv, el is the temperature dependent
!   latent heat of condensation hvap+dldt*(t-ttp), and other values
!   are physical constants defined in parameter statements in the code.
!   The formula is inverted by iterating Newtonian approximations
!   for each the and p until t is found to within 1.e-4 Kelvin.
!   The specific humidity is then computed from pv and pd.
!   This subprogram can be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             exact computation
! 1999-03-01  Iredell             f90 module
!
! Usage:  call stmaxg(tg,the,pk,tma,qma)
!
!   Input argument list:
!     tg         Real(krealfp) guess parcel temperature in Kelvin
!     the        Real(krealfp) equivalent potential temperature in Kelvin
!     pk         Real(krealfp) pressure over 1e5 Pa to the kappa power
!
!   Output argument list:
!     tma        Real(krealfp) parcel temperature in Kelvin
!     qma        Real(krealfp) parcel specific humidity in kg/kg
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp),intent(in):: tg,the,pk
    real(krealfp),intent(out):: tma,qma
    real(krealfp),parameter:: terrm=1.e-4
    real(krealfp) t,p,tr,pv,pd,el,expo,thet,dthet,terr
    integer i
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    t=tg
    p=pk**con_cpor
    do i=1,100
      tr=con_ttp/t
      pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr))
      pd=p-pv
      el=con_hvap+con_dldt*(t-con_ttp)
      expo=el*con_eps*pv/(con_cp*t*pd)
      thet=t*pd**(-con_rocp)*exp(expo)
      dthet=thet/t*(1.+expo*(con_dldt*t/el+el*p/(con_rv*t*pd)))
      terr=(thet-the)/dthet
      t=t-terr
      if(abs(terr).le.terrm) exit
    enddo
    tma=t
    tr=con_ttp/t
    pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr))
    pd=p-pv
    qma=con_eps*pv/(pd+con_eps*pv)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This subroutine computes pressure to the kappa table as a function of pressure
!! for the table lookup function fpkap().
!! Exact pressure to the kappa values are calculated in subprogram fpkapx().
!! The current implementation computes a table with a length
!! of 5501 for pressures ranging up to 110000 Pascals.
  subroutine gpkap
!$$$   Subprogram  documentation  block
!
! Subprogram: gpkap        Compute coefficients for p**kappa
!   Author: Phillips         org: w/NMC2X2   Date: 29 dec 82
!
! Abstract: Computes pressure to the kappa table as a function of pressure
!   for the table lookup function fpkap.
!   Exact pressure to the kappa values are calculated in subprogram fpkapx.
!   The current implementation computes a table with a length
!   of 5501 for pressures ranging up to 110000 Pascals.
!
! Program History Log:
!   94-12-30  Iredell
! 1999-03-01  Iredell             f90 module
! 1999-03-24  Iredell             table lookup
!
! Usage:  call gpkap
!
! Subprograms called:
!   fpkapx     function to compute exact pressure to the kappa
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx
    real(krealfp) xmin,xmax,xinc,x,p
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=0._krealfp
    xmax=110000._krealfp
    xinc=(xmax-xmin)/(nxpkap-1)
    c1xpkap=1.-xmin/xinc
    c2xpkap=1./xinc
    do jx=1,nxpkap
      x=xmin+(jx-1)*xinc
      p=x
      tbpkap(jx)=fpkapx(p)
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This subroutine raises pressure over 1e5 Pa to the kappa power.
!! A linear interpolation is done between values in a lookup table
!! computed in gpkap(). See documentation for fpkapx() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]  p      real, pressure in Pascals
!\param[out] fpkap  real, p over 1e5 Pa to the kappa power
  elemental function fpkap(p)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpkap        raise pressure to the kappa power.
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Raise pressure over 1e5 Pa to the kappa power.
!   A linear interpolation is done between values in a lookup table
!   computed in gpkap. See documentation for fpkapx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy ranges from 9 decimal places
!   at 100000 Pascals to 5 decimal places at 1000 Pascals.
!   On the Cray, fpkap is over 5 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             standardized kappa,
!                                 increased range and accuracy
! 1999-03-01  Iredell             f90 module
! 1999-03-24  Iredell             table lookup
!
! Usage:   pkap=fpkap(p)
!
!   Input argument list:
!     p          Real(krealfp) pressure in Pascals
!
!   Output argument list:
!     fpkap      Real(krealfp) p over 1e5 Pa to the kappa power
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpkap
    real(krealfp),intent(in):: p
    integer jx
    real(krealfp) xj
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xpkap+c2xpkap*p,1._krealfp),real(nxpkap,krealfp))
    jx=min(xj,nxpkap-1._krealfp)
    fpkap=tbpkap(jx)+(xj-jx)*(tbpkap(jx+1)-tbpkap(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function raises pressure over 1e5 Pa to the kappa power.
!! A quadratic interpolation is done between values in a lookup table
!! computed in gpkap(). see documentation for fpkapx() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]   p        real, pressure in Pascals
!\param[out]  fpkapq   real, p over 1e5 Pa to the kappa power
  elemental function fpkapq(p)
!$$$     Subprogram Documentation Block
!
! Subprogram: fpkapq       raise pressure to the kappa power.
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Raise pressure over 1e5 Pa to the kappa power.
!   A quadratic interpolation is done between values in a lookup table
!   computed in gpkap. see documentation for fpkapx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy ranges from 12 decimal places
!   at 100000 Pascals to 7 decimal places at 1000 Pascals.
!   On the Cray, fpkap is over 4 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             standardized kappa,
!                                 increased range and accuracy
! 1999-03-01  Iredell             f90 module
! 1999-03-24  Iredell             table lookup
!
! Usage:   pkap=fpkapq(p)
!
!   Input argument list:
!     p          Real(krealfp) pressure in Pascals
!
!   Output argument list:
!     fpkapq     Real(krealfp) p over 1e5 Pa to the kappa power
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpkapq
    real(krealfp),intent(in):: p
    integer jx
    real(krealfp) xj,dxj,fj1,fj2,fj3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xpkap+c2xpkap*p,1._krealfp),real(nxpkap,krealfp))
    jx=min(max(nint(xj),2),nxpkap-1)
    dxj=xj-jx
    fj1=tbpkap(jx-1)
    fj2=tbpkap(jx)
    fj3=tbpkap(jx+1)
    fpkapq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function raises surface pressure over 1e5 Pa to the kappa power
!! using a rational weighted chebyshev approximation.
!! The numerator is of order 2 and the denominator is of order 4.
!! The pressure range is 40000-110000 Pa and kappa is defined in fpkapx().
!>\param[in]   p        real, surface pressure in Pascals p should be in the 
!! range 40000 to 110000
!\param[out]  fpkapo   real, p over 1e5 Pa to the kappa power
  function fpkapo(p)
!$$$   Subprogram  documentation  block
!
! Subprogram: fpkapo       raise surface pressure to the kappa power.
!   Author: Phillips         org: w/NMC2X2   Date: 29 dec 82
!
! Abstract: Raise surface pressure over 1e5 Pa to the kappa power
!   using a rational weighted chebyshev approximation.
!   The numerator is of order 2 and the denominator is of order 4.
!   The pressure range is 40000-110000 Pa and kappa is defined in fpkapx.
!   The accuracy of this approximation is almost 8 decimal places.
!   On the Cray, fpkap is over 10 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             standardized kappa,
!                                 increased range and accuracy
! 1999-03-01  Iredell             f90 module
!
! Usage:  pkap=fpkapo(p)
!
!   Input argument list:
!     p          Real(krealfp) surface pressure in Pascals
!                p should be in the range 40000 to 110000
!
!   Output argument list:
!     fpkapo     Real(krealfp) p over 1e5 Pa to the kappa power
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpkapo
    real(krealfp),intent(in):: p
    integer,parameter:: nnpk=2,ndpk=4
    real(krealfp):: cnpk(0:nnpk)=(/3.13198449e-1,5.78544829e-2,&
                                         8.35491871e-4/)
    real(krealfp):: cdpk(0:ndpk)=(/1.,8.15968401e-2,5.72839518e-4,&
                                         -4.86959812e-7,5.24459889e-10/)
    integer n
    real(krealfp) pkpa,fnpk,fdpk
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    pkpa=p*1.e-3_krealfp
    fnpk=cnpk(nnpk)
    do n=nnpk-1,0,-1
      fnpk=pkpa*fnpk+cnpk(n)
    enddo
    fdpk=cdpk(ndpk)
    do n=ndpk-1,0,-1
      fdpk=pkpa*fdpk+cdpk(n)
    enddo
    fpkapo=fnpk/fdpk
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function raises pressure over 1e5 Pa to the kappa power.
!! Kappa is equal to rd/cp where rd and cp are physical constants.
!>\param[in]  p        real, pressure in Pascals
!\param[out] fpkapx   real, p over 1e5 Pa to the kappa power
  elemental function fpkapx(p)
!$$$   Subprogram  documentation  block
!
! Subprogram: fpkapx       raise pressure to the kappa power.
!   Author: Phillips         org: w/NMC2X2   Date: 29 dec 82
!
! Abstract: raise pressure over 1e5 Pa to the kappa power.
!   Kappa is equal to rd/cp where rd and cp are physical constants.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   94-12-30  Iredell             made into inlinable function
! 1999-03-01  Iredell             f90 module
!
! Usage:  pkap=fpkapx(p)
!
!   Input argument list:
!     p          Real(krealfp) pressure in Pascals
!
!   Output argument list:
!     fpkapx     Real(krealfp) p over 1e5 Pa to the kappa power
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) fpkapx
    real(krealfp),intent(in):: p
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    fpkapx=(p/1.e5_krealfp)**con_rocp
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This subroutine computes pressure to the 1/kappa table as a function of pressure
!! for the table lookup function frkap().
!! Exact pressure to the 1/kappa values are calculated in subprogram frkapx().
!! The current implementation computes a table with a length
!! of 5501 for pressures ranging up to 110000 Pascals.
  subroutine grkap
!$$$   Subprogram  documentation  block
!
! Subprogram: grkap        Compute coefficients for p**(1/kappa)
!   Author: Phillips         org: w/NMC2X2   Date: 29 dec 82
!
! Abstract: Computes pressure to the 1/kappa table as a function of pressure
!   for the table lookup function frkap.
!   Exact pressure to the 1/kappa values are calculated in subprogram frkapx.
!   The current implementation computes a table with a length
!   of 5501 for pressures ranging up to 110000 Pascals.
!
! Program History Log:
!   94-12-30  Iredell
! 1999-03-01  Iredell             f90 module
! 1999-03-24  Iredell             table lookup
!
! Usage:  call grkap
!
! Subprograms called:
!   frkapx     function to compute exact pressure to the 1/kappa
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx
    real(krealfp) xmin,xmax,xinc,x,p
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=0._krealfp
    xmax=fpkapx(110000._krealfp)
    xinc=(xmax-xmin)/(nxrkap-1)
    c1xrkap=1.-xmin/xinc
    c2xrkap=1./xinc
    do jx=1,nxrkap
      x=xmin+(jx-1)*xinc
      p=x
      tbrkap(jx)=frkapx(p)
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This subroutine raises pressure over 1e5 Pa to the 1/kappa power.
!! A linear interpolation is done between values in a lookup table
!! computed in grkap(). See documentation for frkapx() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]   pkap    real, p over 1e5 Pa to the kappa power
!\param[out]  frkap   real, pressure in Pascals
  elemental function frkap(pkap)
!$$$     Subprogram Documentation Block
!
! Subprogram: frkap        raise pressure to the 1/kappa power.
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Raise pressure over 1e5 Pa to the 1/kappa power.
!   A linear interpolation is done between values in a lookup table
!   computed in grkap. See documentation for frkapx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is better than 7 decimal places.
!   On the IBM, fpkap is about 4 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             standardized kappa,
!                                 increased range and accuracy
! 1999-03-01  Iredell             f90 module
! 1999-03-24  Iredell             table lookup
!
! Usage:   p=frkap(pkap)
!
!   Input argument list:
!     pkap       Real(krealfp) p over 1e5 Pa to the kappa power
!
!   Output argument list:
!     frkap      Real(krealfp) pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) frkap
    real(krealfp),intent(in):: pkap
    integer jx
    real(krealfp) xj
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xrkap+c2xrkap*pkap,1._krealfp),real(nxrkap,krealfp))
    jx=min(xj,nxrkap-1._krealfp)
    frkap=tbrkap(jx)+(xj-jx)*(tbrkap(jx+1)-tbrkap(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function raises pressure over 1e5 Pa to the 1/kappa power.
!! A quadratic interpolation is done between values in a lookup table
!! computed in grkap(). see documentation for frkapx() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]   pkap   real, p over 1e5 Pa to the kappa power
!\param[out]  frkapq real, pressure in Pascals
  elemental function frkapq(pkap)
!$$$     Subprogram Documentation Block
!
! Subprogram: frkapq       raise pressure to the 1/kappa power.
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Raise pressure over 1e5 Pa to the 1/kappa power.
!   A quadratic interpolation is done between values in a lookup table
!   computed in grkap. see documentation for frkapx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is better than 11 decimal places.
!   On the IBM, fpkap is almost 4 times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
!   94-12-30  Iredell             standardized kappa,
!                                 increased range and accuracy
! 1999-03-01  Iredell             f90 module
! 1999-03-24  Iredell             table lookup
!
! Usage:   p=frkapq(pkap)
!
!   Input argument list:
!     pkap       Real(krealfp) p over 1e5 Pa to the kappa power
!
!   Output argument list:
!     frkapq     Real(krealfp) pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) frkapq
    real(krealfp),intent(in):: pkap
    integer jx
    real(krealfp) xj,dxj,fj1,fj2,fj3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xrkap+c2xrkap*pkap,1._krealfp),real(nxrkap,krealfp))
    jx=min(max(nint(xj),2),nxrkap-1)
    dxj=xj-jx
    fj1=tbrkap(jx-1)
    fj2=tbrkap(jx)
    fj3=tbrkap(jx+1)
    frkapq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function raise pressure over 1e5 Pa to the 1/kappa power.
!! Kappa is equal to rd/cp where rd and cp are physical constants.
!>\param[in]   pkap   real, p over 1e5 Pa to the kappa power
!\param[out]  frkapx real, pressure in Pascals
  elemental function frkapx(pkap)
!$$$   Subprogram  documentation  block
!
! Subprogram: frkapx       raise pressure to the 1/kappa power.
!   Author: Phillips         org: w/NMC2X2   Date: 29 dec 82
!
! Abstract: raise pressure over 1e5 Pa to the 1/kappa power.
!   Kappa is equal to rd/cp where rd and cp are physical constants.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   94-12-30  Iredell             made into inlinable function
! 1999-03-01  Iredell             f90 module
!
! Usage:  p=frkapx(pkap)
!
!   Input argument list:
!     pkap       Real(krealfp) p over 1e5 Pa to the kappa power
!
!   Output argument list:
!     frkapx     Real(krealfp) pressure in Pascals
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) frkapx
    real(krealfp),intent(in):: pkap
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    frkapx=pkap**(1/con_rocp)*1.e5_krealfp
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This subroutine computes lifting condensation level temperature table
!! as a function of temperature and dewpoint depression for function ftlcl().
!! Lifting condensation level temperature is calculated in subprogram ftlclx()
!! The current implementation computes a table with a first dimension
!! of 151 for temperatures ranging from 180.0 to 330.0 Kelvin
!! and a second dimension of 61 for dewpoint depression ranging from
!! 0 to 60 Kelvin.
  subroutine gtlcl
!$$$     Subprogram Documentation Block
!
! Subprogram: gtlcl        Compute equivalent potential temperature table
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute lifting condensation level temperature table
!   as a function of temperature and dewpoint depression for function ftlcl.
!   Lifting condensation level temperature is calculated in subprogram ftlclx
!   The current implementation computes a table with a first dimension
!   of 151 for temperatures ranging from 180.0 to 330.0 Kelvin
!   and a second dimension of 61 for dewpoint depression ranging from
!   0 to 60 Kelvin.
!
! Program History Log:
! 1999-03-01  Iredell             f90 module
!
! Usage:  call gtlcl
!
! Subprograms called:
!   (ftlclx)    inlinable function to compute LCL temperature
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    integer jx,jy
    real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,tdpd,t
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xmin=180._krealfp
    xmax=330._krealfp
    ymin=0._krealfp
    ymax=60._krealfp
    xinc=(xmax-xmin)/(nxtlcl-1)
    c1xtlcl=1.-xmin/xinc
    c2xtlcl=1./xinc
    yinc=(ymax-ymin)/(nytlcl-1)
    c1ytlcl=1.-ymin/yinc
    c2ytlcl=1./yinc
    do jy=1,nytlcl
      y=ymin+(jy-1)*yinc
      tdpd=y
      do jx=1,nxtlcl
        x=xmin+(jx-1)*xinc
        t=x
        tbtlcl(jx,jy)=ftlclx(t,tdpd)
      enddo
    enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
!> This function computes temperature at the lifting condensation level
!! from temperature and dewpoint depression.
!! A bilinear interpolation is done between values in a lookup table
!! computed in gtlcl(). See documentation for ftlclx() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]  t      real, LCL temperature in Kelvin
!>\param[in]  tdpd   real, dewpoint depression in Kelvin
!\param[out] ftlcl  real, temperature at the LCL in Kelvin
  elemental function ftlcl(t,tdpd)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftlcl        Compute LCL temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute temperature at the lifting condensation level
!   from temperature and dewpoint depression.
!   A bilinear interpolation is done between values in a lookup table
!   computed in gtlcl. See documentation for ftlclx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is better than 0.0005 Kelvin.
!   On the Cray, ftlcl is ? times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
! 1999-03-01  Iredell             f90 module
!
! Usage:   tlcl=ftlcl(t,tdpd)
!
!   Input argument list:
!     t          Real(krealfp) LCL temperature in Kelvin
!     tdpd       Real(krealfp) dewpoint depression in Kelvin
!
!   Output argument list:
!     ftlcl      Real(krealfp) temperature at the LCL in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftlcl
    real(krealfp),intent(in):: t,tdpd
    integer jx,jy
    real(krealfp) xj,yj,ftx1,ftx2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xtlcl+c2xtlcl*t,1._krealfp),real(nxtlcl,krealfp))
    yj=min(max(c1ytlcl+c2ytlcl*tdpd,1._krealfp),real(nytlcl,krealfp))
    jx=min(xj,nxtlcl-1._krealfp)
    jy=min(yj,nytlcl-1._krealfp)
    ftx1=tbtlcl(jx,jy)+(xj-jx)*(tbtlcl(jx+1,jy)-tbtlcl(jx,jy))
    ftx2=tbtlcl(jx,jy+1)+(xj-jx)*(tbtlcl(jx+1,jy+1)-tbtlcl(jx,jy+1))
    ftlcl=ftx1+(yj-jy)*(ftx2-ftx1)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function computes temperature at the lifting condensation level
!! from temperature and dewpoint depression.
!! A biquadratic interpolation is done between values in a lookup table
!! computed in gtlcl(). see documentation for ftlclx() for details.
!! Input values outside table range are reset to table extrema.
!>\param[in]   t      real, LCL temperature in Kelvin
!>\param[in]   tdpd   real, dowpoint depression in Kelvin
!\param[out]  ftlcl  real, temperature at the LCL in Kelvin
  elemental function ftlclq(t,tdpd)
!$$$     Subprogram Documentation Block
!
! Subprogram: ftlclq       Compute LCL temperature
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute temperature at the lifting condensation level
!   from temperature and dewpoint depression.
!   A biquadratic interpolation is done between values in a lookup table
!   computed in gtlcl. see documentation for ftlclx for details.
!   Input values outside table range are reset to table extrema.
!   The interpolation accuracy is better than 0.000003 Kelvin.
!   On the Cray, ftlclq is ? times faster than exact calculation.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
! 1999-03-01  Iredell             f90 module
!
! Usage:   tlcl=ftlclq(t,tdpd)
!
!   Input argument list:
!     t          Real(krealfp) LCL temperature in Kelvin
!     tdpd       Real(krealfp) dewpoint depression in Kelvin
!
!   Output argument list:
!     ftlcl      Real(krealfp) temperature at the LCL in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftlclq
    real(krealfp),intent(in):: t,tdpd
    integer jx,jy
    real(krealfp) xj,yj,dxj,dyj
    real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33
    real(krealfp) ftx1,ftx2,ftx3
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    xj=min(max(c1xtlcl+c2xtlcl*t,1._krealfp),real(nxtlcl,krealfp))
    yj=min(max(c1ytlcl+c2ytlcl*tdpd,1._krealfp),real(nytlcl,krealfp))
    jx=min(max(nint(xj),2),nxtlcl-1)
    jy=min(max(nint(yj),2),nytlcl-1)
    dxj=xj-jx
    dyj=yj-jy
    ft11=tbtlcl(jx-1,jy-1)
    ft12=tbtlcl(jx-1,jy)
    ft13=tbtlcl(jx-1,jy+1)
    ft21=tbtlcl(jx,jy-1)
    ft22=tbtlcl(jx,jy)
    ft23=tbtlcl(jx,jy+1)
    ft31=tbtlcl(jx+1,jy-1)
    ft32=tbtlcl(jx+1,jy)
    ft33=tbtlcl(jx+1,jy+1)
    ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21
    ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22
    ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23
    ftlclq=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function computes temperature at the lifting condensation level
!! from temperature and dewpoint depression.  the formula used is
!! a polynomial taken from Phillips mstadb routine which empirically
!! approximates the original exact implicit relationship.
!>\param[in]    t    real, temperature in Kelvin
!>\param[in]    tdpd real, dewpoint depression in Kelvin
!\param[out]   ftlclo   real, temperature at the LCL in Kelvin
  function ftlclo(t,tdpd)
!$$$   Subprogram  documentation  block
!
! Subprogram: ftlclo       Compute LCL temperature.
!   Author: Phillips         org: w/NMC2X2   Date: 29 dec 82
!
! Abstract: Compute temperature at the lifting condensation level
!   from temperature and dewpoint depression.  the formula used is
!   a polynomial taken from Phillips mstadb routine which empirically
!   approximates the original exact implicit relationship.
!   (This kind of approximation is customary (inman, 1969), but
!   the original source for this particular one is not yet known. -MI)
!   Its accuracy is about 0.03 Kelvin for a dewpoint depression of 30.
!   This function should be expanded inline in the calling routine.
!
! Program History Log:
!   91-05-07  Iredell             made into inlinable function
! 1999-03-01  Iredell             f90 module
!
! Usage:  tlcl=ftlclo(t,tdpd)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!     tdpd       Real(krealfp) dewpoint depression in Kelvin
!
!   Output argument list:
!     ftlclo     Real(krealfp) temperature at the LCL in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftlclo
    real(krealfp),intent(in):: t,tdpd
    real(krealfp),parameter:: clcl1= 0.954442e+0,clcl2= 0.967772e-3,&
                                    clcl3=-0.710321e-3,clcl4=-0.270742e-5
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ftlclo=t-tdpd*(clcl1+clcl2*t+tdpd*(clcl3+clcl4*t))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This function computes temperature at the lifting condensation level
!! from temperature and dewpoint depression.  A parcel lifted
!! adiabatically becomes saturated at the lifting condensation level.
!! The water model assumes a perfect gas, constant specific heats
!! for gas and liquid, and neglects the volume of the liquid.
!! The model does account for the variation of the latent heat
!! of condensation with temperature.  The ice option is not included.
!! The Clausius-Clapeyron equation is integrated from the triple point
!! to get the formulas:
!!\n  pvlcl=con_psat*(trlcl**xa)*exp(xb*(1.-trlcl))
!!\n  pvdew=con_psat*(trdew**xa)*exp(xb*(1.-trdew))
!!\n where pvlcl is the saturated parcel vapor pressure at the LCL,
!! pvdew is the unsaturated parcel vapor pressure initially,
!! trlcl is ttp/tlcl and trdew is ttp/tdew.  The adiabatic lifting
!! of the parcel is represented by the following formula
!!\n    pvdew=pvlcl*(t/tlcl)**(1/kappa)
!!\n This formula is inverted by iterating Newtonian approximations
!! until tlcl is found to within 1.e-6 Kelvin.  Note that the minimum
!!   returned temperature is 180 Kelvin.
!>\param[in]  t      real, temperature in Kelvin
!>\param[in]  tdpd   real, dewpoint depression in Kelvin
!\param[out] ftlclx real, temperature at the LCL in Kelvin
  elemental function ftlclx(t,tdpd)
!$$$   Subprogram  documentation  block
!
! Subprogram: ftlclx       Compute LCL temperature.
!   Author: Iredell          org: w/NMC2X2   Date: 25 March 1999
!
! Abstract: Compute temperature at the lifting condensation level
!   from temperature and dewpoint depression.  A parcel lifted
!   adiabatically becomes saturated at the lifting condensation level.
!   The water model assumes a perfect gas, constant specific heats
!   for gas and liquid, and neglects the volume of the liquid.
!   The model does account for the variation of the latent heat
!   of condensation with temperature.  The ice option is not included.
!   The Clausius-Clapeyron equation is integrated from the triple point
!   to get the formulas
!       pvlcl=con_psat*(trlcl**xa)*exp(xb*(1.-trlcl))
!       pvdew=con_psat*(trdew**xa)*exp(xb*(1.-trdew))
!   where pvlcl is the saturated parcel vapor pressure at the LCL,
!   pvdew is the unsaturated parcel vapor pressure initially,
!   trlcl is ttp/tlcl and trdew is ttp/tdew.  The adiabatic lifting
!   of the parcel is represented by the following formula
!       pvdew=pvlcl*(t/tlcl)**(1/kappa)
!   This formula is inverted by iterating Newtonian approximations
!   until tlcl is found to within 1.e-6 Kelvin.  Note that the minimum
!   returned temperature is 180 Kelvin.
!
! Program History Log:
! 1999-03-25  Iredell
!
! Usage:  tlcl=ftlclx(t,tdpd)
!
!   Input argument list:
!     t          Real(krealfp) temperature in Kelvin
!     tdpd       Real(krealfp) dewpoint depression in Kelvin
!
!   Output argument list:
!     ftlclx     Real(krealfp) temperature at the LCL in Kelvin
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
    real(krealfp) ftlclx
    real(krealfp),intent(in):: t,tdpd
    real(krealfp),parameter:: terrm=1.e-4,tlmin=180.,tlminx=tlmin-5.
    real(krealfp) tr,pvdew,tlcl,ta,pvlcl,el,dpvlcl,terr,terrp
    integer i
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    tr=con_ttp/(t-tdpd)
    pvdew=con_psat*(tr**con_xpona)*exp(con_xponb*(1.-tr))
    tlcl=t-tdpd
    do i=1,100
      tr=con_ttp/tlcl
      ta=t/tlcl
      pvlcl=con_psat*(tr**con_xpona)*exp(con_xponb*(1.-tr))*ta**(1/con_rocp)
      el=con_hvap+con_dldt*(tlcl-con_ttp)
      dpvlcl=(el/(con_rv*t**2)+1/(con_rocp*tlcl))*pvlcl
      terr=(pvlcl-pvdew)/dpvlcl
      tlcl=tlcl-terr
      if(abs(terr).le.terrm.or.tlcl.lt.tlminx) exit
    enddo
    ftlclx=max(tlcl,tlmin)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end function
!-------------------------------------------------------------------------------
!> This subroutine computes all physics function tables.  Lookup tables are
!! set up for computing saturation vapor pressure, dewpoint temperature,
!! equivalent potential temperature, moist adiabatic temperature and humidity,
!! pressure to the kappa, and lifting condensation level temperature.
  subroutine gfuncphys
!$$$     Subprogram Documentation Block
!
! Subprogram: gfuncphys    Compute all physics function tables
!   Author: N Phillips            w/NMC2X2   Date: 30 dec 82
!
! Abstract: Compute all physics function tables.  Lookup tables are
!   set up for computing saturation vapor pressure, dewpoint temperature,
!   equivalent potential temperature, moist adiabatic temperature and humidity,
!   pressure to the kappa, and lifting condensation level temperature.
!
! Program History Log:
! 1999-03-01  Iredell             f90 module
!
! Usage:  call gfuncphys
!
! Subprograms called:
!   gpvsl       compute saturation vapor pressure over liquid table
!   gpvsi       compute saturation vapor pressure over ice table
!   gpvs        compute saturation vapor pressure table
!   gtdpl       compute dewpoint temperature over liquid table
!   gtdpi       compute dewpoint temperature over ice table
!   gtdp        compute dewpoint temperature table
!   gthe        compute equivalent potential temperature table
!   gtma        compute moist adiabat tables
!   gpkap       compute pressure to the kappa table
!   grkap       compute pressure to the 1/kappa table
!   gtlcl       compute LCL temperature table
!
! Attributes:
!   Language: Fortran 90.
!
!$$$
    implicit none
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    call gpvsl
    call gpvsi
    call gpvs
    call gtdpl
    call gtdpi
    call gtdp
    call gthe
    call gtma
    call gpkap
    call grkap
    call gtlcl
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  end subroutine
!-------------------------------------------------------------------------------
end module
